Curso 6 Inferencia estadística

VOLVER

Curso 6 Inferencia estadística

INFERENCIA ESTADÍSTICA
TEORÍA DE LA PROBABILIDAD
DISTRIBUCIÓN DE PROBABILIDAD
PROBABILIDAD CONDICIONAL
ESTUDIO DETALLADO DE LA SENSIBILIDAD Y ESPECIFICIDAD
VARIABLES INDEPENDIENTES
VALOR ESPERADO O ESPERANZA MATEMATICA(Centro geométrico de las Áreas de probabilidad) o MEDIA
VARIABILIDAD
MEDIA Y DESVIACIÓN ESTANDAR (ERROR ESTANDAR) DE MÚLTIPLES MUESTRAS
EJEMPLOS CON DISTRIBUCIONES
ESTUDIO ASINTÓTICO
INTERVALOS DE CONFIANZA
CONTRASTE DE HIPÓTESIS (Hypothesis testing)
P-VALUE
TEST
PODER DE CONTRASTE
COMPARACIONES MÚLTIPLES (MULTIPLE TESTING)
RESAMPLING INFERENCE

INFERENCIA ESTADÍSTICA

La inferencia estadística es el conjunto de técnicas estadísticas que nos permiten llegar a conclusiones sobre una población a partir de los datos de una muestra de dicha población. La estadística por tanto parte de los datos para extraer unas conclusiones, es un método deductivo. Para dar soporte teórico a la estadística existe el cálculo de probabilidades, que lo que proporciona es un método inductivo, es decir desde una regla es capaz de llegar a la misma conclusión. Tanto para la estadística como para la probabilidad, existen dos corrientes:

Frecuentista: Se basa en la repetición continuada de los hechos. Cuanto más se repita un experimento, la probabilidad de que ocurra un suceso será más regular y su frecuencia relativa se acercará mucho a la probabilidad del suceso.

USAREMOS ESTE MÉTODO CUANDO NO SEPAMOS ANTES LAS PROBABILIDADES DE CADA SUCESO, POR EJEMPLO QUE UNA MONEDAÉ ESTÉ TRUCADA Y NO SEPAMOS CÓMO, SÓLO POREMOS SABER SU PARCIALIDAD REPITIENDO, POR ELLO SE LLAMA CONOCIMIENTO A POSTERIORI

La interpretación clásica, mayoritaria por lo menos hasta ahora, define la probabilidad en términos de experimentación. Si repites un experimento un número infinito de veces y compruebas que en 350 de cada 1.000 ocasiones se ha producido un determinado resultado, un frecuentista diría que la probabilidad de ese resultado es del 35%. Basándose en esta definición, un frecuentista afirma que es posible asociar a cada evento una probabilidad de obtener un valor VERDADERO del mismo.

La aproximación clásica se basa por lo tanto en estudiar la probabilidad real de las cosas, tratando de determinar hasta qué punto una medición realizada sobre un conjunto de experimentos se aproxima a la probabilidad real que subyace. Es por ello que un frecuentista definiría probabilidad como una expresión matemática que predice en qué medida es verosímil que ciertos eventos ocurran basándose en el patrón observado hasta este momento.

El enfoque frecuentista considera que los datos son aleatorios debido a la naturaleza del proceso que los genera y que cada vez que repitamos el experimento obtendremos un resultado diferente. De la misma forma cree que la hipótesis será cierta o falsa para el proceso estudiado pero que debido a la naturaleza aleatoria de los datos nuestro análisis puede señalarla como falsa cuando es cierta (falso negativo) o verdadera cuando es falsa (falso positivo). Por ello el enfoque frecuentista nos proporcionará, además de una respuesta a la pregunta de si la hipótesis es cierta o falsa, una probabilidad denominada p-valor. Este número es un indicador cuantas veces tendríamos que repetir la toma de datos y el análisis estadístico para obtener el resultado opuesto. Nótese que, al ser los datos aleatorios, siempre puede ser que acabe usando un conjunto que me apunte en la dirección contraria a la que realmente se comporta la población total.

Bayesiana: Hace un anális basado en la certeza que se tiene de que un hecho ocurrirá. Las observaciones se usan para actualizar o inferir la posibilidad de una hipótesis. Es una estrategia complementaria a la frecuentista, ya que con cada iteración añadiríamos un conocimiento obtenido en esa iteración para enriquecer el cálculo de la probabilidad.

USAREMOS ESTE MÉTODO CUANDO SEPAMOS A PRIORI LAS PROBABILIDADES DE CADA SUCESO, POR EJEMPLO AL LANZAR UNA MONEDA PERFECTA, POR ESO SE LLAMA CONOCIMIENTO A PRIORI

Por el contrario, la interpretación bayesiana se basa en un conocimiento limitado de las cosas. Afirma que sólo asocias una probabilidad a un evento porque hay incertidumbre sobre el mismo, es decir, porque no conoces todos los hechos. En realidad, un evento dado, o bien ocurrirá (probabilidad=100%) o bien no ocurrirá (probabilidad=0%). Cualquier otra cosa es una aproximación que hacemos del problema a partir de nuestro conocimiento incompleto del mismo. El enfoque bayesiano se basa por lo tanto en la idea de refinar predicciones a partir de nuevas evidencias. Un bayesiano definiría probabilidad como la expresión matemática que mide el nivel de conocimiento que tenemos para hacer una predicción. Por lo tanto, para un bayesiano, estríctamente hablando es incorrecto decir predigo que hay un 30% de probabilidades de que ocurra el evento P, sino que debería decir basándome en el conocimiento actual que tengo, tengo un 30% de certeza de que P ocurrirá.

Los frecuentistas piensan en los parámetros como unos valores fijos pero desconocidos. La estimación se basa en la elección de aquellos valores de los parámetros que maximizan la probabilidad de observar los datos.

Por su parte, los bayesianos interpretan los parámetros como variables aleatorias cuya distribución de probabilidad es estudiada en base al Teorema de Bayes. La idea es simple, un bayesiano ha de tener una distribución subjetiva de los parámetros antes de ver los datos (a priori) que modificará en función de los datos que haya observado para obtener una distribución a posteriori que resumirá todo el conocimiento del investigador sobre los parámetros de interés dados los datos y sus creencias a priori.

Los bayesianos le dan la vuelta al enfoque frecuentista. Para ellos los datos son fijos, no aleatorios porque ¿no se han obtenido ya? ¿Qué tienen de aleatorios los valores anotados en esas tablas que hemos ido generando? En cambio la hipótesis es para ellos aleatoria y puede ser o no verdad. Los bayesianos no determinan si la hipótesis es cierta o falsa, únicamente calculan la probabilidad (en sentido bayesiano) de que sea cierta o falsa.

Finalmente, sabemos que hay muchos resultados frecuentistas que pueden obtenerse desde una perspectiva bayesiana aunque la filosofía subyacente sea diferente. Mínimos cuadrados ordinarios (MCO) es el ejemplo más sencillo. Dicho estimador frecuentista coincide exactamente con la media de la distribución bayesiana bajo unas creencias concretas en el marco del modelo lineal normal. Existen en la literatura muchos otros ejemplos.

El objetivo de las técnicas de inferencia estadística es contrastar la hipótesis nula (Ho) y la hipótesis alternativa (H1), y en este sentido se pueden cometer dos tipos de error:

Error tipo I o alfa --> Rechazar la hipótesis nula siendo esta verdadera.
Error tipo II o beta --> Rechazar la hipótesis alternativa siendo esta cierta. La probabilidad complementaria del tipo II se llama poder de contraste (p-value) y es 1-beta.
                                     Este valor está totalmente impactado por el tamaño de la muestra, ya que el intervalo de confianza se va reduciendo cuanto mayor es el espacio muestral, por tanto más probable es rechazar Ho de manera efectiva.

la base de la inferencia estadistica es la probabilidad. Esta, asigna un número entre 0 y 1 que significa la posibilidad de que sea cierto. La probabilidad es el modelo de inferencia que se usa habitualmente para hacer estimaciones sobre fenómenos aparentemente aleatorios.

Teniendo en cuenta la conexión entre ambos enfoques en la práctica y la posibilidad de obtener los mismos resultados bajo ambos enfoques, los bayesianos argumentan que su interpretación del problema es siempre más intuitiva y natural. Recordemos que un bayesiano proporcionaría conclusiones del tipo: hay un 95% de probabilidad de que el parámetro esté entre 0,3 y 0,8, sin embargo, un frecuentista afirmaría: si generáramos 100 muestras aleatorias del mismo tamaño y repitiéramos la estimación 100 veces, en 95 de ellas el parámetro estimado estaría entre 0,3 y 0,8.

Los contrastes de hipótesis parecen por tanto más naturales en el marco bayesiano. Para un bayesiano convencido, un p-valor o un t-ratio es una herejía, ya que lo único que necesitamos para contrastar hipótesis es tener la distribución a posteriori de los parámetros. Salvo en contadas excepciones, un frecuentista basa su contraste de hipótesis en el análisis asintótico, es decir, en calcular p-valores de la distribución asintótica del estimador (no del parámetro verdadero, que es un valor fijo). Esta distribución, mayoritariamente normal gracias a numerosos teoremas centrales del límite, es la que el estimador tendría si el investigador tuviera muchos más datos de los que realmente tiene. ¿es ésta una buena forma de hacer inferencia en general? Por el contrario, el bayesiano se basa en la distribución de los parámetros dados sus datos, que puede en principio tener cualquier forma no gaussiana.

Al admitir información adicional, los métodos bayesianos reconocen que cada problema es distinto y promueven que el procedimiento de análisis se adapte al problema en cuestión, y no al revés; como consecuencia, tienden a ser más flexibles. Sin embargo, para poder hacer esto posible, la implementación de las técnicas bayesianas usualmente requiere de un esfuerzo computacional muy alto. La mayor parte de este esfuerzo se concentra en el cálculo de ciertas características de la distribución final del parámetro de interés. Por ejemplo, hay que integrar para pasar de una distribución conjunta a una colección de distribuciones marginales que sean útiles para hacer inferencias sobre los parámetros de interés. En la gran mayoría de los problemas las integrales requeridas no pueden resolverse analíticamente, por lo que es necesario contar con métodos numéricos eficientes que permitan calcular o aproximar integrales en varias dimensiones.

El enfoque bayesiano no sustituye el estudio de las frecuencias, lo enriquece con experiencias anteriores ya que puede ir refinando hacia atrás los resultados obtenidos en cada iteración. Un ejemplo extremo es que si en los últimos 4.500 millones de años en Cuenca ha salido el sol por el levante, suponemos que la probabilidad de que mañana vaya a volver a salir el sol es más alta con un estudio bayesiano de lo que sería con un análisis de estadística clásica, que no da ningún peso a la historia pasada. Pero cuando un médico pronostica le quedan unos seis meses de vida, querríamos creer que el galeno nos está facilitando esta información no sólo en base a un análisis de las frecuencias para esta dolencia, sino que además incorpora su experiencia de otros diagnósticos similares, y en especial nuestro propio historial clínico. Las conclusiones son afectadas por la premisa, pero si ésta es bien formulada, más que un defecto sería una ventaja

Afortunadamente, durante la segunda mitad del siglo XX se comenzaron a desarrollar técnicas numéricas flexibles y eficientes basadas en métodos de simulación estocástica. Esto, aunado al desarrollo de la tecnología que ha dado lugar a una mayor capacidad de procesamiento y de almacenamiento de los equipos de cómputo gracias a la virtualización, el Cloud Computing y el Big Data, está provocando un auge de los métodos bayesianos.

Probablidad es la cantidad de población que resume una aleatoriedad.Dos herramientas muy útiles para el cálculo de probabilidad son la distribución de probabilidad y la función de distribución

Esto es, ante la pregunta de si una hipótesis es correcta, los frecuentistas responden un sí o un no pero apostillan con qué frecuencia obtendrías el resultado opuesto mientras que los bayesianos no responden directamente y se limitan a indicar cómo de seguros están de uno u otro resultado.

Haremos usos de ambos métodos dependiendo de lo que necesitemos calcular. Usaremos por ejemplos métodos frecuentistas para el cálculo del error muestral, mientra que usaremos métodos bayesianos para el cálculo de especificidad y sensibilidad.

TEORÍA DE LA PROBABILIDAD

qownnotes-media-FbSrHM

Se asigna una probabidad dependiendo del conjunto (set) en el que se quiera colocar

qownnotes-media-uHCHMM

La suma de todas las posibilidades excuyentes es 1

qownnotes-media-yNdmps

qownnotes-media-yWCfyx

qownnotes-media-GBLpWb

qownnotes-media-waZsrq

qownnotes-media-VafxWO

qownnotes-media-mTAmEg

\[ \begin{eqnarray*} A_1 & = & \{\mbox{Person has sleep apnea}\} \\ A_2 & = & \{\mbox{Person has RLS}\} \end{eqnarray*} \] \[ \begin{eqnarray*} P(A_1 \cup A_2 ) & = & P(A_1) + P(A_2) - P(A_1 \cap A_2) \\ & = & 0.13 - \mbox{Probability of having both} \end{eqnarray*} \]

qownnotes-media-cIgzcE

Un ejemplo de variable discreta es todo lo que sea contar un número de cosas(por ejemplo machos y hembras, o cualquier otra variable categórica), un ejemplo de continua es la medida de tu peso durante un mes.

Supongamos que deseamos conocer algún parámetro poblacional, por ejemplo la edad media de una población. En el contexto de la estadística clásica (frecuentista), podríamos extraer al azar múltiples y diferentes muestras representativas de forma sucesiva y determinar el estimador más apropiado en cada una de esas muestras. Para este caso, el estimador más apropiado para la media poblacional es la media muestral. Cada una de las medias de estas muestras tomarían un determinado valor diferente para cada muestra (o lo que es lo mismo, el estimador es tratado como una variable en sí mismo).En estadística, un estimador es un estadístico (esto es, una función de la muestra) usado para estimar un parámetro desconocido de la población.

Suponiendo la normalidad de la variable edad, los valores que irían tomando estas medias muestrales se distribuirían según una curva de Gauss (algunos de estos valores se darían con mayor frecuencia y otros con menor frecuencia). El área bajo dicha curva de distribución de frecuencias del estimador es la probabilidad total, es decir, la unidad.

DISTRIBUCIÓN DE PROBABILIDAD

-Función de probabilidad (función de masa de probabilidad, PMF en ingles) Función que determina la probabilidad de cada evento posible, cuando los posibles valores son DISCRETOS: ejemplo:

valor p(1) -> 0.3 valor p(4) ->0.6 valor p(7)->0.1

PMF valores discretos

Las distribuciones de probabilidad en base a su PMF más comunes son Bernoulli (valor p o p-1 muy útil para la prevalencia de una enfermedad), binomial (bernuilli repetido n veces) o poisson(contar el número de veces que se da un evento en el tiempo)

-Función de densidad (PDF en ingles) Función que determina la probabilidad de cada evento posible, cuando los posibles valores son CONTINUOS, por ejemplo probabilidad asociada a tiempo de espera en urgencias.

Ejemplo la proporción de llamadas que pueden ser atendidas en un día sigue la siguiente curva:

\[ f(x) = \left\{\begin{array}{ll} 2 x & \mbox{ for } 1 > x > 0 \\ 0 & \mbox{ otherwise} \end{array} \right. \]

x <- c(-0.5, 0, 1, 1, 1.5); y <- c( 0, 0, 2, 0, 0)
plot(x, y, lwd = 3, frame = FALSE, type = "l")

Cuál será la probabilidad de atender el 75% al menos de las llamadas?

Esta operación es muy sencilla ya que se puede calcular con trigonometría o através de la función beta, que es justo para distribuciones de probabiidad trianangular

1.5 * .75 / 2

## [1] 0.5625

pbeta(.75, 2, 1)

## [1] 0.5625

Por tanto la probabilidad de atender al menos el 75% de las llamadas es del 56%

Ya que calculamos el área bajo la curva, si elegimos un valor concreto su pdf nos dara 0,ya que el area de una recta es 0.

Area bajo curva de una recta es 0

Las distribuciones de probabilidad en base a su PDF más comunes son Normal (la normal estandar es la que tiene e(x) = 0 y d(x) = 1)o de Gaus o exponencial (el tiempo entre dos eventos)

Los histogramas, pueden representar esta función de densidad en base a probabilidades o bien en base a número de ocurrencias de cada evento. EN nuestro caso, esto es fundamental ya que las probabilidades que vamos a manejar en el mundo de la la sanidad áestn todas basadas en la ocurrencia de un suceso, no en un cálculo a priori de sus probabilidades.

http://www.statmethods.net/graphs/density.html

FUNCION DE DISTRIBUCION(CDF en ingles) Es la probabilidad igual o menor asociada a un valor (probabilidad acumulada). Nos va a resultar muy útil para hacer el cálculo de la media en variables contínuas, ya que para su cálculo necesitamos este valor asociado a la variable.

\[ F(x) = P(X \leq x) \] Por otro lado la función de supervivencia será la que represente la probabilidad mayor a un valor

\[ S(x) = P(X > x) \]

Como se puede ver por tanto $S(x) = 1 - F(x)$

En variables contnuas, la PDF es la derivada de la CDF, como se puede ver en el ejemplo:

\[ F(x) = P(X \leq x) = \frac{1}{2} Base \times Height = \frac{1}{2} (x) \times (2 x) = x^2 \] La función de supervivencia por tanto será: \[ S(x) = 1 - x^2 \]

Esto puede ser calculado facilmente con la función pbeta, por ejmplo la probabilidad de atender el 40,50 y 60% de las llamadas

pbeta(c(0.4, 0.5, 0.6), 2, 1)

## [1] 0.16 0.25 0.36

Las probabilidades de la CDF se puede expresar también en Cuantiles y Percentiles

CUANTILES cuando tienes cuantil 0.95 significa que el 95% de la población está por debajo de tí, es decir se centra en el valor de la probabilidad menor o igual en lugar del resultado.

El cuantil alfa es el valor de la X que te da un valor concreto de cuantil. Ejemplo Cuantil alfa de 0.98 es X=6 PERCENTIL es el cuantil expresado en porcentaje

La media es el percentil 50%

En este ejemplo, queremos saber el cuantil 0.5 es decir la probabilidad de atender el 50% de las llamadas $0.5 = F(x) = x^2$

sqrt(0.5)

## [1] 0.7071068

Podemos hacer uso téambin de nuevo de pbeta:

qbeta(0.5, 2, 1)

## [1] 0.7071068

Hasta aqui hemos hablado en todo momento de medias sin tener en cuenta los datos, esto es así porque estamos calculando la media de la POBLACIÓN, un modelo de probabilidad conecta los datos con la población a través de asunciones. Por tanto la media de la que hemos estado hablando es el estimando, mientras que la media de la máuestra será el estimador.

PROBABILIDAD CONDICIONAL

La probabilidad de que un suceso pase sabiendo que ha pasado otro.

qownnotes-media-UsnVoJ

Un ejemplo de probabilidad condicionada, es qué probabilidades tengo de sacar un 1 si sí que ha salido impar.

qownnotes-media-jvqfcq

Teorema de Bayes Es el uso más famoso de la probabilidad condicional. Calcula la probabilidad de A condicionado a B desde la probabilidad de que pase B condicionado a A.

qownnotes-media-bfGADK

Es muy útil para calcular la efectividad de test diagnósticos. Ejemplo: + y - son los resultados de diagnóstico de un test médico, y D que el sujeto tenga la enfermedad.

SENSIBILIDAD es la probabilidad de que el test sea positivo sabiendo que individuo tiene la enfermedad. Es desable que sea alto.

qownnotes-media-sYMoUf

ESPECIFIDAD es la probabilidad de que el test sea negativo sabiendo que el individuo no tiene la enfermedad. Tambi?n es desable que sea alto

Prevalencia aparente–> Probabilidad de tener sepsis observada. Prevalencia real –> es practicamente imposible de determinar pero se puede aproximar desde la prevalencia aparente con el estimador Rogan-Gladen

Sensibilidad -> Proporci?n de casos acertados con alerta positiva

Especificidad -> Proporci?n de casos acertados con alerta negativa

Valor predictivo positivo –> Probabilidad de que si me hago un test sea positivo (tiene en cuenta prevalencia)

Valor predictivo negativo –> Probabilidad de que si me hago un test sea negativo (tiene en cuenta prevalencia)

Razón de verosimilitud -> Existen con el fin de dimensionar el beneficio clínico de un test de manera independiente a la prevalencia, ya que esta puede variar. El ratio entre la posibilidad de observar un paciente séptico frente a un paciente no séptico

$DLR_+ = .997 / (1 - .985) \approx 66$ LR+ Ratio de 124 significa que la probabilidad de tener la enfermedad después del test es 66 veces mayor que antes del test

$DLR_- = (1 - .997) / .985 \approx .003$ LR- es 0.005 significa que la probabilidad de no tener la enfermedad antes del test es 0.003 veces mayor que después del test

Odds ratio -> Asociación entre exposición y resultado. Ejemplo odds ratio de indice de cancer de pulmón en fumador es 2. Esto significa que fumadores tienen 2 veces más probab. de tener cancer.

Diferencia entre probability y odds:

Probabilidad 0.80 Odds 0.80/1-0.80 = 4 –> Odds de 4 a 1

Test	Séptico	No séptico	TOTAL
Pos	8503	2758	11261
Neg	37	355024	355061
TOTAL	8540	357782	366322

Prevalencia aparente= (8503+37)/(8503+37+2758+355024)=0.023

Prevalencia real = (Prev.apa+Espec-1)/(Sensib.+Espec-1) = (0.023+0.992-1)/(0.995+0.992-1)= 0.01519757

Proporción falso positivo = 2758/11261 = 0.24

Proporción falson negativo = 37/355061 = 0.0001

Sensibilidad = 8503/(8503+37) = 0.995

Especificidad = 355024/(2758+355024) = 0.992

VPP = 8503/(8503+2758) = 0.75

VPN = 355024/(2578+355024) = 0.992

LR+ = (sensibilidad/1-especificidad) = 124

LR- = (1-sensibilidad/especificidad) = 0,005

Odds ratio (raz?n de oportunidades) = (8503/2758)/(37/355024) = 29582.43

Ratio o razón –> Indice obtenido al dividir dos cantidades. Ninguno de los elementos del numerador estón incluidos en el denominador

Ejemplo Ratio alerta posit = 8503/2758 = 3.08 por cada 3 alertas                            salta 1 que no es

Proporción -> El numerador est? incluido en el denominador

    Proporcion alerta posit =  8503/(8503+2758) = 0.75

Tasa -> Magnitud de cambio de un parámetro por unidad de cambio de otro que incluye una medida adiconal de tiempo en el denominador

Tasa de pacientes septicos en 2018: 8540/366322 = 0.023 es una tasa
de 23 pacientes septicos por cada 1000 urgencias.

qownnotes-media-tBAWRb

qownnotes-media-cYnYCe

POSITIVE PREDICTIVE VALUE es la probabilidad de tener la enfermedad sabiedo que el test es positivo.

NEGATIVE PREDICTIVE VALUE es la probabilidad de no tener la enfermedad sabiendo que el test es negativo.

En la ausencia de un test, se llama PREVALENCIA a la probabilidad de tener esa enfermedad.

qownnotes-media-dOUpkP

qownnotes-media-tLsDNl

qownnotes-media-pPYGcB

qownnotes-media-RUIsxz

DLR+ es cada cuantos test positivos acertados se da un test negativo err?neo en este caso 998/15 = 66, cada 66 test positivos acertados se da 1 negativo err?neo

DRL- es cada cuantos test positivos erróneos se da un test negativo acertado.

qownnotes-media-PaNiRC

qownnotes-media-xTKrgI

qownnotes-media-CtTfPs

qownnotes-media-DynpEw

qownnotes-media-htWhsk

qownnotes-media-EtEyvi

ESTUDIO DETALLADO DE LA SENSIBILIDAD Y ESPECIFICIDAD

Una vez que la prueba se ha aplicado a una población determinada, varios índices se utilizan para evaluar la prueba: - Verdadero positivos (VP): Número de casos que la prueba declara positivos y que son verdaderamente positivos. - Falsos positivos (FP): Número de casos que la prueba declara positivos y que en realidad son negativos. - Verdadero negativo (VN): Número de casos que la prueba declara negativos y que son realmente negativos. - Falsos negativos (FN): Número de casos que la prueba declara negativos y que en realidad son positivos. - N = TP + FP + FN + TN

qownnotes-media-HqjvcR

Sensibilidad (equivalente a la tasa de positivos verdaderos): Proporci?n de casos positivos que están bien detectadas por la prueba. La definición matemática es: Sensibilidad = VP / (VP + FN).
Especificidad (también llamada Tasa de verdaderos negativos): proporci?n de casos negativos que son bien detectadas por la prueba. La definición matemática es: Especificidad = VN / (VN + FP).
Tasa de falsos positivos (FPR): Proporci?n de casos negativos que la prueba detecta como positivos.
Falso Negativo Precio (FNR): Proporci?n de casos positivos que la prueba detecta como negativo. - Prevalencia: la frecuencia relativa de los acontecimiento de inter?s en la muestra total (VP + FN) / N.
Valor Predictivo Positivo (VPP): Proporci?n de casos verdaderamente positivos entre los casos positivos detectados por la prueba. Tenemos PPV = TP / (TP + FP), o PPV = Sensibilidad x Prevalencia / [(Sensibilidad x Prevalencia + (1-Especificidad) (1-Prevalencia)]. Es un valor fundamental que depende de la prevalencia, un índice que es independiente de la calidad de la prueba.
Valor predictivo negativo (VPN): Proporción de casos verdaderamente negativos entre los casos negativos detectados por la prueba. Tenemos VPN = VN / (VN + FN), o VPP = Especificidad x (1 - Prevalencia) / [(Especificidad (1-Prevalencia) + (1-sensibilidad) x Prevalencia] . Este ?ndice depende tambi?n de la prevalencia que se independiente de la calidad de la prueba.
Razón de verosimilitud positiva (LR +): Esta relación indica a qué punto una persona tiene más posibilidades de ser positivo en la realidad cuando la prueba es que est? diciendo es positivo. Tenemos LR + = sensibilidad / (1-especificidad). La RP + es un valor positivo o nulo.
Razón de verosimilitud Negativa (LR-): Esta relación indica a qué punto una persona tiene más posibilidades de ser negativo, en realidad, cuando la prueba es que est? diciendo es positivo. Tenemos LR-= (1-sensibilidad) / (especificidad). El LR-es un valor positivo o nulo.
Odds ratio: El odds ratio indica la cantidad de un individuo es m?s probable que sea positivo si el resultado es positivo, en comparaci?n con los casos en que la prueba es negativa. Por ejemplo, un odds ratio de 2 significa que la probabilidad de que el caso positivo se produce es dos veces superior si la prueba es positiva que si es negativo. El odds ratio es un valor positivo o nulo. Tenemos Odds ratio = VPxVN / (FPxFN).
El riesgo relativo: El riesgo relativo es un ratio que mide mejor la prueba se comporta cuando se trata de un informe positivo que cuando es negativo. Por ejemplo, un riesgo relativo de 2 significa que la prueba es dos veces más potente cuando es positiva que cuando es negativo. Un valor cercano a 1 corresponde a un caso de independencia entre las filas y columnas, y una prueba de que funciona tan bien cuando es positiva que cuando es negativo. El riesgo relativo es un valor nulo o un valor positivo dado por: riesgo relativo = VP / (VP + FP) / (FN / (FN + VN)).

VARIABLES INDEPENDIENTES

Dos variables son independientes si la ocurrencia de una no tiene nada que ver en la ocurrencia de otra.

Por ejemplo tirar dos veces una moneda y que salga las dos cara, son sucesos independientes 0.5x0.5 = 0.25

Si los sucesos no son independientes la intersección de las probabilidades no se pueden calcular multiplicando.

qownnotes-media-edrncZ

Todos los estudios que realizamos no sólo están basados en sucesos independientes, si no que además están identicamente distribuidos, IID variables

qownnotes-media-GvDchr

VALOR ESPERADO O ESPERANZA MATEMATICA(Centro geométrico de las Áreas de probabilidad) o MEDIA

Se insiste mucho en la diferencia entre muestra y población. En este caso cuando hagamos cálculos sobre una muestra haremos los cálculos en base a los valores encontrados en esa muesta (Ej:tiro un dado 10 veces y salen 2 1s, 3 2s, 1 3, 4 5s.). Cuando se habla sobre los cáculos sobre una población estaremos hablando sobre las probabilidades que tenga una variable de tomar un valor sin haber realizado muestras ( o lo que es lo mismo si hubiéramos hecho infinitas muestras)

MEDIA DE UNA MUESTRA Es la suma de todos los valores dividido por el número de muestras TirarDado = {1,1,5,4,3,4,6,5,4,6} _ x = (1+1+5+4+3+4+6+5+4+6)/10 = 3.9

MEDIA DE UNA POBLACIÓN Es una propiedad de las distribuciones de probabilidad de la misma manera que lo son la mediana (valor del punto que está en la mitad) y la varianza

En variables discretas será el valor de su probabilidad multiplicado por el propio valor

E(x)= x1P(x1)+x2P(x2)+…

$$
E[X] = 1 \times \frac{1}{6} + 2 \times \frac{1}{6} +
3 \times \frac{1}{6} + 4 \times \frac{1}{6} +
5 \times \frac{1}{6} + 6 \times \frac{1}{6} = 3.5
$$

qownnotes-media-S12388

En variables contínuas, tenemos que hacer uso de la función de distribución para recoger la probabilidad acumulada.(dx)

E(x) = Integral (xf(x) d(x))

En PDF es lo mismo , el valor esperado es el centro de masa del área que encierra su función de densidad. Una propieda de la media es que su valor es el mismo en la función de densidad (DE POBLACION) que en la función de densidad de la media de X muestras del mismo tipo (MUESTRA) , es decir, si el valor esperado de una distribución normal centrada en 0 es 0 (en azul), función de densidad de la media de 10 distribuciones normales distintas será también 0 (en rosa). Cuanto mayor sea la población sobre la que estamos sacando la media (en lugar de 10 que fuera 100) más se aglutinarán las densidades entorno al valor esperado.

qownnotes-media-z14800

MEDIA DE MÚLTIPLES MUESTRAS

La media de un conjunto de muestras es en sí también una variable aleatoria, de la que podríamos calcular su media y su varianza ( la media de las medias y la varianza de las medias). _ E(x) = suma x / n

VARIABILIDAD

VARIANZA DE UNA MUESTRA Es la suma de la distancia de todos los puntos de la muestra respecto a la media (se pone al cuadrado para evitar distancias negativas y calcular su distancia absoluta)

S^2 = (1-3.9 +1-3.9 +5-3.9 +4-3.9 +3-3.9 +4-3.9 +6-3.9 +5-3.9 +4-3.9 +6-3.9 )^2 / 10 -1 = 0

S = Desviación estándar –> Se aplica a raiz cuadrada para que esté en las mismas unidades que el resto.

VARIANZA DE UNA POBLACIóN Varianza: es la dispersión de los datos respecto a la media.

Var(X) = E(X-mu)^2 (Mide la distancia de todos los puntos con respecto a la media, al cuadrado para que sean distancias absolutas) sd(X) = sqrt(var(x)) Desviación estándar es la raíz cuadrada de la varianza de manera que el dato está en las mismas unidades que el eje de las X, lo cual es muy útil.

desviación típica = desviación estándar = error estándar = error típico

También tiene una distribución de población asociada así como un valor esperado llamado varianza de la población. La desviación estandar se llama también error estandar.

En el caso del ejemplo:Var(X) = E(x)^2 - mu^2

E(x)^2 = (11/6)^{2+(21/6)^{2+(31/6)^2+(41/6)}2+(51/6)}2+(61/6)^2 = 15.17

mu^2 = 3.5^2 = 12,25

Var(X) = 15.17-12.25 = 2.92

qownnotes-media-dfmfMN

En el caso de la varianza de las muestras hay una pequeña variación ya que se divide por n-1. Esto es debido a que las primeras n-1 observaciones son son independientes de la media. la úiltima no es independiente ya que puede ser calculada por medio de la media de la muestra usada en la formula. En otras palabras, la varianza de la muestra es CASI la media de la desviación al cuadrado de la media de la muestra. _ Var(x) suma var(x) /n-1

La distribución de la varianza de las muestras también se verá afectado por el tamaño de la muestra, de manera que a mayor tamaño de la muestra la curva se irá acercando a la variaza de la población que está intentando estimar:

qownnotes-media-M19368

Var(X’)=Var(1/n*Sum(X_i))=(1/n^{2)Var(Sum(X_i))=(1/n^2)Sum(sigma}2)=sigma^2/n para poblaciones infinitas cuando las muestras son independientes | infinite populations when our samples are independe

VARIANZA DE MÚLTIPLES MUESTRAS

La desviación estándar de una una estadística se llama error estándar y es la varianza dividido el número de muestras, o lo que es lo mismo la desviación estándar dividido raíz cuadrada del número de muestas (numerador y denominador aplicados la raíz cuadrada)

desviación estándar = S^2/n = S/raizCuadrada(n)

qownnotes-media-CiZnlf

qownnotes-media-A14800

MEDIA Y DESVIACIÓN ESTANDAR (ERROR ESTANDAR) DE MÚLTIPLES MUESTRAS

Si simulamos con R 1000 muestras de una distribción normal con 10 valores posibles vemos que se acerca mucho a lo calculado para la población, y cuanto más se aumenten las muestras más cerca estará ese valor:

qownnotes-media-RlgzCm

qownnotes-media-hYQmGV

qownnotes-media-KWuTdE

qownnotes-media-MkxXYP Hay una medida adicional que es el coeficiente de dispersión de Pearson, que es la desviación típica dividida por su media, de esta manera podemos medir la dispersión de dos variables de distintas unidades. Por ejemplo, comparar la dispersóin del peso de los niños de una clase frente a la altura de otra clase.

El índice de Gini también nos puede servir para medir cómo de concentrados están los datos

El coeficiente de Fisher nos dice cómo de asimétrica es la curva de distribución.

La courtosis es cómo de aglutinados están las frecuencias alrededor de la media.

EJEMPLOS CON DISTRIBUCIONES

BERNULLI

Si X es una variable aleatoria que mide el número de éxitos, y se realiza un único experimento con dos posibles resultados (éxito o fracaso), se dice que la variable aleatoria X se distribuye como una Bernoulli de parámetro p. (Ejemplo tirar 1 moneda 1 vez y que salga cruz)

qownnotes-media-MnhxZU

Si hacemos varias observaciones de Bernuilli pongamos $x_1,\ldots, x_n$ , la posibilidad (likelihood)de que el evento sea éxito es \[ \prod_{i=1}^n p^{x_i} (1 - p)^{1 - x_i} = p^{\sum x_i} (1 - p)^{n - \sum x_i} \] $\hat p = \sum_i x_i / n$ es el máximo estimador de posibilidad para $p$

Todas las posibilidades para n pequeña

BINOMIAL

cuenta el número de éxitos en una secuencia de n ensayos de Bernoulli independientes entre sí, con una probabilidad fija p de ocurrencia del éxito entre los ensayos. Un experimento de Bernoulli se caracteriza por ser dicotómico, esto es, solo dos resultados son posibles

Para n = 1, la binomial se convierte, de hecho, en una distribución de Bernoulli.

qownnotes-media-zTNfQz

qownnotes-media-dnOuDU

qownnotes-media-mjITSs

qownnotes-media-fjsoAv

plot(pvals, dbinom(7, 8, pvals) / dbinom(7, 8, 7/8) ,
     lwd = 3, frame = FALSE, type = "l", xlab = "p", ylab = "likelihood")

DISTRIBUCION NORMAL ESTANDAR

Es una distribución normal pero con la particularidad de que tiene mu = 0 y sigma = 1

qownnotes-media-BjqEKS

De manera que 1 desviación estandar recoge el 68% de los datos, 2 desviaciones estándar el 95% y 3 desviaciones estandar el 99%.

qownnotes-media-SIFoSu

qownnotes-media-hjbGgE

DISTRIBUCION NORMAL

Por tanto una distribucion normal tendrá mu != 0 y sigma != 1, la única diferencia con los gráficos anteriores es que no estará centrada en 0, y que los ejes de las x cada 1,2 y 3 estar? multiplicado por su sigma.

qownnotes-media-mAfrpi

qownnotes-media-kOrKhD

qownnotes-media-RxbLpG

qownnotes-media-R14800

qownnotes-media-B14800

qownnotes-media-P19368

qownnotes-media-L19368

POISSON

expresa, a partir de una frecuencia de ocurrencia media, la probabilidad de que ocurra un determinado n?mero de eventos durante cierto periodo de tiempo. Concretamente, se especializa en la probabilidad de ocurrencia de sucesos con probabilidades muy peque?as, o sucesos “raros”.

Tanto el valor esperado como la varianza de una variable aleatoria con distribución de Poisson son iguales

qownnotes-media-mMWyDG

qownnotes-media-BmoCZj

qownnotes-media-Y14800

qownnotes-media-i14800

qownnotes-media-i19368

qownnotes-media-u19368

Es muy útil para:

*Conteos de datos Cuantas veces se repite un determinado evento.

*Modelar eventos de tiempo o de supervivencia, que modeliza el tiempo que se tarda en que ocurra un determinado suceso. Por el nombre de la técnica parecería que se analizara el tiempo hasta la muerte (Análisis de supervivencia) pero, en realidad, puede analizarse cualquier otro suceso.

En Análisis de supervivencia la muestra consiste en el seguimiento de una serie de individuos desde el inicio del estudio hasta su final y, ante una situación de este tipo, es frecuente que se produzca la desaparición de alguno de esos individuos que entran en el estudio. También es posible que al entrar un individuo en el estudio, éste termine antes de que en ese individuo se produzca el suceso que se pretende detectar. Un ejemplo es el tiempo en que cominezan a aparecer determinados sintomas en pacientes sanos: qownnotes-media-gtvKdN

*Para modelar tablas de contingencia.Modelando tablas de contingencia (cruzar dos características y ver el número de casos que pertenecen a cada categoria)

qownnotes-media-uSdaBj

*Aproximaciones binomiales cuando n es muy grande y p muy pequeño, por ejemplo en epidemiología el tener o no tener una enfermedad sobre una población.

qownnotes-media-p14800

qownnotes-media-l14800

qownnotes-media-n19368

qownnotes-media-m19368

ESTUDIO ASINTÓTICO

Es el término para el análisis del comportamiento de las estadísticas cuando el tamaño de la muestra tiene un límite (normalmente infinito). Es la base de estudio frecuentista de probabilidades.

Un estimador es consistente si converge con lo que queremos estimar. La ley de números largos (LLN) dice que la media de la muestra, la varianza y la desviación estandar de iid muestras es consistente con la media de la población, la varianza y la desviación estándar respectivamente.

Teorema central del limite (CLT) : dice que la distribución de las medias siguen una distribución normal estándar centrada en el origen (ya que estamos restando la media de la población), con una desviación estándar igual al error estándar de la media, pero debe ser para un número de repeticiones suficiente como para que se cumpla.

$X_1,\ldots,X_n$ es una colección de variables aleatorias iid con media $\mu$ y varianza $\sigma^2$ $\bar X_n$ es la media de la muestra Entonces $\frac{\bar X_n - \mu}{\sigma / \sqrt{n}}$ tiene una distribución normal para largos $n$.

Por ejemplo, si tiramos un dado muchas veces y hacemos n observaciones

$\frac{\bar X_n - 3.5}{1.71/\sqrt{n}}$

Se ve claramente la dierencia de hacer 1,2 y 6 tiradas 10000 veces. La curva normal estándar lógicamente va a estar centrada en el origen ya que le estamos restando la media.

Al seguir una distribución normal, podemos usar los cuantiles para saber con qué probabilidad la media obtenida estará incluida en ese intervalo, para una distrióbucin normal sabemos \[\bar X_n \pm z_{1-\alpha/2}\sigma / \sqrt{n}\] y lo llamaremos intervalo de confianza.

En este ejemplo de las alturas de los hijos, podemos dar un intervalo de confianza en el cuál con un 95% de posibilidades la media obtenida en la muestra está contenida en la media de la población:

library(UsingR);data(father.son); x <- father.son$sheight

## Loading required package: MASS

## Loading required package: HistData

## Loading required package: Hmisc

## Loading required package: lattice

## Loading required package: survival

## Loading required package: Formula

## Loading required package: ggplot2

##
## Attaching package: 'Hmisc'

## The following objects are masked from 'package:base':
##
##     format.pval, units

##
## Attaching package: 'UsingR'

## The following object is masked from 'package:survival':
##
##     cancer

(mean(x) + c(-1, 1) * qnorm(.975) * sd(x) / sqrt(length(x))) / 12

## [1] 5.709670 5.737674

Esta misma operación la podemos hacer con Bernuilli para proporciones de éxito o fallo:

\[ \hat p \pm z_{1 - \alpha/2} \sqrt{\frac{p(1 - p)}{n}} \] Sabiendo que la moneda es imparcial, p = 1/2, p(1-p)=1/4 y sqrt(p(1-p)) = 1/2

Aplicando le intervalor Wald del 95%

$p(1-p) \leq 1/4$ for $0 \leq p \leq 1$ $\alpha = .05$ so that $z_{1 -\alpha/2} = 1.96 \approx 2$ \[ 2 \sqrt{\frac{p(1 - p)}{n}} \leq 2 \sqrt{\frac{1}{4n}} = \frac{1}{\sqrt{n}} \] Entonces $\hat p \pm \frac{1}{\sqrt{n}}$ es un buen intervalo de confianza rápido para $p$

Un ejemplo, es que en un sondeo de 100 personas, 56 votan sí, tengo la votación ganada? 1/sqrt(100)=.1 esto nos da un intervalo de confianza de Wald (95%) de entre (0.46, 0.66) así que no tenemos nada claro.

(0.56 + c(-1, 1) * qnorm(.975) * sqrt(0.56*0.44/100))

## [1] 0.4627099 0.6572901

binom.test(56,100)$conf.int

## [1] 0.4571875 0.6591640
## attr(,"conf.level")
## [1] 0.95

Si el proceso no fuera parcial (p=0.9 por ejemplo) el CTL nos garantiza únicamente que es una buena aproximación en el límite, pero se puede observar que si hacemos la misma prueba que con p = 0.5 tardan muchas más iteraciones en converger hacia una normal estándar.

Hay veces que el intervalo de wald no se comporta de manera correcta si el número de repeticiones no es el adecuado, y para evitarlo se puede aplicar la corrección de Agresti/Coull añadiendo dos aciertos y dos errores $p = X_n +2 / n+4 $

Por último si se trata de una distribución de Poisson por ejemplo para conteos:

Vemos que damos un un error diagnóstico 5 veces cada 94.32 días, cuál es el porcentaje de fallo diario con un intervalo de confianza dentro del 95%?

$X \sim Poisson(\lambda t)$ Estimate $\hat \lambda = X/t$ $Var(\hat \lambda) = \lambda / t$

\[ \frac{\hat \lambda - \lambda}{\sqrt{\hat \lambda / t}} = \frac{X - t \lambda}{\sqrt{X}} \rightarrow N(0,1) \]

x <- 5
t <- 94.32
lambda <- x/t
round(lambda + c(-1, 1) * qnorm(.975) * sqrt(lambda/t),3 )

## [1] 0.007 0.099

poisson.test(5,94.32)$conf.int

## [1] 0.01721254 0.12371005
## attr(,"conf.level")
## [1] 0.95

En este caso podemos observar que el intervalo calculado es más conservador ya que t no es lo suficientemente grande.

qownnotes-media-bG7876

Hay veces que el número de repeticiones no es suficiente en distribuciones discretas binomiales y de Poisson y se puede hacer una aproximación del valor estimado sumando +2 al numerador y +4 al denominador (añadiendo 2 exitos y 2 fallos)

qownnotes-media-jY7876

qownnotes-media-ZUuUAB

INTERVALOS DE CONFIANZA

Resumiendolo que hemos visto anteriormente, la media de los valores medios se aproxima a una normal estandar, con media $\mu$ y sd $}{\sigma / \sqrt{n}}$.

Por tanto, al seguir una normal estandar, el desplazamiento desde la media de 2 unidades del error estándar nos dará un intervalo en el que tendremos la probabilidad del 95% de contener el valor medio. A este intervalo lo llamamos ZQuantile X SE

$\bar X_n \pm 2\sigma / \sqrt{n}$ es el llamado intervalo 95% de $\mu$

En la vida real, las muestras con las que trabajamos no suelen ser tan grandes, por lo que necesitamos una adaptación, utilizando la distribución T Student de Gosset. surge del problema de estimar la media de una población normalmente distribuida cuando el tamaño de la muestra es pequeño. Está también centrada en 0. Si no fuera simétrica en 0 (skewed distr) habría que calculalo de otra manera o coger otro estimador como la mediana

qownnotes-media-WYnisx

Aparece de manera natural al realizar la prueba t de Student para la determinación de las diferencias entre dos medias muestrales y para la construcción del intervalo de confianza para la diferencia entre las medias de dos poblaciones cuando se desconoce la desviación típica de una población y ésta debe ser estimada a partir de los datos de una muestra.

qownnotes-media-MIqXgZ

En adelante vamos a ver pequeñas muestras de datos y en este caso tendremos que usar la distribución T de Gosset y intervalos de confianza T, la única diferencia con la ecuación anterior es que usaremos el TQuantile

Cuando una distribución T creece y crece en muestras se acaba convirtiendo en distribución Z.

La distribución T tiene los extremos más pequños que la normal.A diferencia de la dist estandar, que tiene 2 parametros (media y varianza). La T de gosset tiene sólo un parámetro llamado grados de libertad t = n-1, cuanto más crezca este número más se parecerá a una standar normal, ya que si n es muy grande la diferencia entre una estandar normal, que es dividida por sigma, y la T de student que es dividida por la varianza, se hace insignificativa, es algo similar a lo que nos pasaba al estudiar la varianza de las muestras, que cuanto mayor era el número de las muestras más se acerca al teórico de la población. En este caso para calcular T con n-1 grados de libertad:

qownnotes-media-KN1524

El valor de esta t nos dará un intervalo de confianza, para el cuál dependiendo del número de individuos, en 95% de los casos ese valor está incluido, en el caso de que quisiéramos usar ese quantil qt(0.975,n)

Dependiendo de los grados de confianza, la diferencia de las distribuciones será mayor, en este caso los grados de confianza son muy amplios, lo que hace que sea muy similar

qownnotes-media-ASAszo

La distribución T es muy útil cuando queremos comparar cómo varía un grupo frente a otro, pero siempre teniendo en cuenta que deben ser datos iid normales, centrados en su media

Vamos a coger dos grupos y vamos a estudiar la diferencia de variabilidad entre ellos, primero estando relaccionados (que los individuos sean los mismos, y queremos ver su evolución) y luego siendo grupos independientes (Son grupos de personas totalmente distintas) .

Esa diferencia de variabilidad estará encerrada en un intervalo de confianza y se suele utilizar el intervalo T para compararlos

Los resultados del estudio por tanto varían dependiendo de si se emparejan los resultados porque están relaccionados o son independientes: Los datos son un grupo de 10 pacientes a los que se les ha tratado con 2 drogas distintas y se intenta comparar los resultados

data(sleep)
sleep

##    extra group ID
## 1    0.7     1  1
## 2   -1.6     1  2
## 3   -0.2     1  3
## 4   -1.2     1  4
## 5   -0.1     1  5
## 6    3.4     1  6
## 7    3.7     1  7
## 8    0.8     1  8
## 9    0.0     1  9
## 10   2.0     1 10
## 11   1.9     2  1
## 12   0.8     2  2
## 13   1.1     2  3
## 14   0.1     2  4
## 15  -0.1     2  5
## 16   4.4     2  6
## 17   5.5     2  7
## 18   1.6     2  8
## 19   4.6     2  9
## 20   3.4     2 10

En las siguientes gráficas se ven cada uno de los valores de ambos resultados, estando la leyenda la correspondencia entre el color de la raya que une a los individuos y el número de individuo que es en el grupo

DEPENDIENTES (paired=TRUE)

## [1] -3.363874  0.203874

qownnotes-media-JY1524

qownnotes-media-rV1524

Si queremos calcular el intevalo de confianza de los resultados emparejados, creamos un vector nuevo con la diferencia de los resultados de ambos grupos difference <- g2-g1

y se calcula la media y la desviación estandar mn <- mean(difference) s <- sd(difference) Los calculos se pueden hacer a mano mn + c(1,-1)qt(.975)(s/sqrt(10)) o con la función t.test()

qownnotes-media-qs1524

qownnotes-media-ay1524

Esto lo que nos dice, es que con 95% de posibilidades ( la funcion qt(.975) nos da la t para el 95% con 9 grados de libertad) la diferencia media de los efectos entre ambas drogas es entre 0.7 y 2.46 horas de sueño

data(sleep)
x1 <- sleep$extra[sleep$group == 1]
x2 <- sleep$extra[sleep$group == 2]
n1 <- length(x1)
n2 <- length(x2)
sp <- sqrt( ((n1 - 1) * sd(x1)^2 + (n2-1) * sd(x2)^2) / (n1 + n2-2))
md <- mean(x1) - mean(x2)
semd <- sp * sqrt(1 / n1 + 1/n2)
md + c(-1, 1) * qt(.975, n1 + n2 - 2) * semd

## [1] -3.363874  0.203874

t.test(x1, x2, paired = FALSE, var.equal = TRUE)$conf

## [1] -3.363874  0.203874
## attr(,"conf.level")
## [1] 0.95

t.test(x1, x2, paired = TRUE)$conf

## [1] -2.4598858 -0.7001142
## attr(,"conf.level")
## [1] 0.95

INDEPENDIENTES (PAIRED= FALSE)

En este caso tenemos dos grupos distintos de personas, unos que tomaron el medicamento y otro que tomaron placebo y vamos a medir su presión arterial. No se pueden emparejar porque son personas distintas y los grupos pueden tener incluso tamaño diferente.

qownnotes-media-Fz1524

El objetivo de nuevo es encontrar el intervalo de confianza de la diferencia media de la presión arterial media entre ambos grupos con un 95% de posibilidad. Para hacer el cálculo de los intervalos de confianza vamos a trabajar restando la media de ambos conjuntos, el quantile se calcula en base a ambos conjuntos, ya que hemos calculado su media por separado, así que hay que restar 2 a la suma de ambos para obtener el grado de confianza de todo.

qownnotes-media-aV1524

La varianza común es más dificil de calcular, haremos una estimación media entre ambas varianzas aplicando pesos en cada una de ellas dependiento del tamaño de cada grupo diviendolo por el tamaño total. Se ve más claro en este ejemplo extraído del libro Fundamentals of Biostatistics. Primero haremos los cálculos suponiendo que las varianzas son equivalentes. Se puede observar que el intervalo de confianza sobre difrencia media de ambas presiones arteriales medias contiene el valor 0, por lo que no se puede descartar que las medias de ambos grupos sean iguales.

qownnotes-media-Tn1524

Aquí se puede apreciar el error de cálculo si tratamos el ejemplo anterior como grupos independientes

qownnotes-media-kP1524

Se puede apreciar en el cálculo de ambos casos que cuando los tratamos de manera distinta un cálculo contiene el 0 y el otro no lo contiene, lo que significa que cuando los grupos están emparejados tienen muchas más medias distintas, de manera opuesta a cuando no los emparejamos y los tratamos como independientes en que se contiene el valor 0 y por tanto las diferencias en las medias son menos distintas ya que se tratan todos los datos como un solo grupo.

La función t.test tiene otro parámetro con el que se puede jugar que es si la varianza de los datos es igual en ambos conjuntos. En este ejemplo, se intenta hacer un intervalo de confianza entre la ganancia de peso de pollos con 4 dietas, comparándolas de 2 en 2, en este ejmplo la 1 y la 4, cuya diferencia de variablidad es muy evidente a nivel gráfico.

qownnotes-media-Bt1524

qownnotes-media-mA1524

vemos que los resultados de las estimaciones varían dependiendo del tratamiento de la relacción entre las varianzas de los grupos.

qownnotes-media-ll1524

DISTINTA VARIANZA ENTRE GRUPOS

En el caso de que la varianza de los grupos no sea la misma o no estemos seguros, En caso de duda, se debe utilizar varianza no es igual, para calcularala

qownnotes-media-Mp1524

qownnotes-media-Rh1524

CHI SQUARE TEST

Este test se utiliza para saber si dos variables son independientes, es decir, si la aparición de una está relaccionada (ojo no es la causa, sólo si están relacconada) con la otra. Es apropiado utilizar este método en las siguientes situaciones:

El muestreo es aleatorio
las variables son categoricas
si los datos recogidos en la tabla de contingencia tienen en las frecuencias esperadas al menos 5 casos.

Este método consiste en 4 pasos:

Establecer la hipótesis: Ho-> Variable A y variable B son independientes Ha-> Variable A y variable B no son independientes La hipótesis alternativa por tanto es que sabiendo el valor A nos puede ayudar a predecir el valor de B Formular el plan de análisis Este paso describe cómo usar muestras de los datos para aceptar o rechazar la hipótesis nula: Nivel de significancia elegir un nivel apropiado entre 0 y 1 Test a utilizar en este caso chi-cuadrado Analizar la muestra de los datos Grados de libertad DF = (r-1)(c-1) donde r son el nº de niveles de A y c el nº de niveles de B Frecuencia esperada = Son computadas por separado para cada nivel de cada variable en cada nivel de la otra variable Er,c = (nr*nc)/n Donde Er,c es la frecuencia esperada para el nivel r de la variable A y nivel c de la B nr es el número total de observaciones del nivel r en A nc es el número total de observaciones del nivel c en B n es el número total de la muestra Test estadístico = x2 = sum[(Or,c - Er,c)2/Er,c] Donde Or,c es la frecuencia del nivel r de la variable A y nivel C de la B Er,c es la frecuencia esperada del nivel r de la variable A y el nivel c de la B P-value = es la probabilidad de observar una muestra de manera aleatoria al menos tan extrema como la nuestra. Para poder conocer esta probabilidad, se debe calcular si la distribución pertenece a una chi-cuadrado, si es así, utilizar los grados de libertad anteriores y el test estadístico para obtener el pvalue de la muestra de datos.

CONTRASTE DE HIPÓTESIS (Hypothesis testing)

Un contraste de hipótesis (también denominado test de hipótesis o prueba de significación) es un procedimiento para juzgar si una propiedad que se supone en una población estadística es compatible con lo observado en una muestra de dicha población.La aplicación de cálculos probabilísticos permite determinar a partir de qué valor debemos rechazar la hipótesis garantizando que la probabilidad de cometer un error es un valor conocido a priori.

A partir de una muestra de la población en estudio, se extrae un estadístico (media muestral, varianza muestral,etc.) cuya distribución de probabilidad esté relacionada con la hipótesis en estudio y sea conocida. Se toma entonces como región de rechazo al conjunto de valores que es más improbable bajo la hipótesis, esto es, el conjunto de valores para el que rechazaremos la hipótesis nula si el valor del estadístico observado entra dentro de él.

Para probar una hipótesis vamos a comparar siempre con la hipótesis nula Ho, e invalidarla, ya que partimos de que la hipótesis nula es siempre cierta.

A continuación, vamos a ver cuatro ejemplos de cómo descartar H_0, con una distribución normal grande, con una distribución normal pequeña ambos con varianza equvalente, uno con varianza no equivalente y un ejemplo con distribición binomial.

Vamos a trabajar con un conjunto de 100 pacientes con sobrepeso que tienen desorden del sueño y una media de 32 eventos/hora durante la noche con una desviación estándar de 10 eventos/hora. Se considera que desorden del sueño es 30 eventos/hora. Vamos a utilizar el Test-Z como test estadístico, para poder usarlo necesitamos la media , la varianza y la distribución.

Queremos invalidar H0 que dice que mu = 30 y validar H_a que dice que mu > 30 (32 en este caso concreto)

ERROR TIPO I –> EL investigador SI encuentra diferencia entre Ho y Ha cuando no la hay (Falso positivo) ALFA probabilidad de cometer ese error

ERROR TIPO II –> El investigador NO encuentra diferencia entre Ho y Ha cuando sí la hay (falso negativo) BETA probabilidad de cometer ese error

Error tipo I –> Rechazar Ho cuando es correcta esto es alfa (probabilidad de cometer este error es típicamente 5% llamada alfa) Tampoco queremos un alfa muy bajo ya que nunca rechazaríamos H_o aunque fuera falsa. Error tipo II –> Aceptar H_o cuando es falso Error tipo III –> Rechazar H_a cuando es correcta Error tipo IV –> Aceptar H_a cuando es falsa

qownnotes-media-amRzlL

El error estándar de la media de un conjunto es desviacion_estandar/sqrt(población), en este ejemplo 10/sqrt(100) = 1

H_o por tanto tiene un media normalmente distribuida,mu = 30 y varianza = 1 (estimada como el error estandar que es 1 al cuadrado)

Estamos buscando una constante C que hiciera que la probabilidad de que la media sea mayor que C es %5 Esto es P(X>C|H_o) = 5%

qownnotes-media-k12128

El área sin colorear se puede calcular como qnorm(.95)

Recordando que en una distribución normal el percentil 95 está a 1.645 desviaciones estándar de la media y teniendo una distribución normal de media 30 y 1 de desviación estandar N(30,1)

Rechazaríamos la hipotesis nula si la media fuera mayor que un valor (C) cuya probabilidad sea 5% C = 30 + 1x1.645 = 31.645

qownnotes-media-Kv7296

Rechazar la hipótesis nula cuando la media es >= 31.645 tiene una probabilidad de fallo del 5%, nuestra media es 32 que es mayor y cae dentro de ese área así que rechazamos H_o En este gráfico vemos las dos distribuciones, y cómo la media de H_a está dentro del área del 5% de H_o, lo que nos sirve para rechazar que nuestro resultado sea por azar en un 95% de los casos al menos.

qownnotes-media-T16176

El Z-score es lo mismo pero calculando cuantos errores estándar la media de la hipótesis es diferente de la media, esto se calcula computando la distancia entre las medias (32-30) y dividiendo por el error estandar de la media (s/sqrt(n)) que en este caso es 10/sqrt(100)

Z-score = 2, significa que está a 2 errores estandar de la media.

Un alumno A saca una puntuación de 85 en un examen cuyas puntuaciones tienen una media de 79 con una desviación típica de 8. Un alumno B saca 74 en un examen cuyas puntuaciones tienen una media de 70 y desviación típica de 5. ¿Cuál de los dos alumnos obtuvo una puntuación mejor? La respuesta, desde el punto de vista de la unidad tipificada, se obtiene así:

Las puntuaciones tipificadas de los alumnos A y B son respectivamente:

qownnotes-media-blJJQy

Así el alumno B lo hizo mejor que el A, aunque su puntuación de 74 es inferior a 85.

qownnotes-media-rn7296

La regla por tanto para descartar H_0 es (X-mu)*sqrt(n)/s > Z

Digamos que ahora bajamos de 100 ejemplos a 16, se convierte entonces en una distribución T con 15 grados de libertad. Aquí ya no podemos utlizar Zscore como test estadístico porque la distribución es normal, pero con pocos individuos (distribución T-student)

qownnotes-media-Hy7296

En este caso fallamos al rechazar Ho porque 0.8< 1.7531

En lugar de concentrar todo el error en uno de los extremos, se suele distribuir a ambos lados

qownnotes-media-sv7296

Normalmente todos estos cálculos no los tenemos que hacer a mano, ya que la función t.test() nos da toda esta información

qownnotes-media-nS7296

En este ejemplo, queremos saber si la media de altura de los hijos es equivalente a la media de altura de los padres. Cogemos los datos de las alturas de padres e hijos, y restamos la altura de un padre a la de su hijo, y vamos a ver si la diferencia (haciendo la resta se asume que están emparejados) es 0 o no es 0 con t.test (tambien se podrían pasar ambos arrays a la función y como tercer parametro paired.equals=TRUE)

Ambos t(cuantas desviaciones estandar está la media del conjunto del de H_o) y df (grados de libertad) son muy grandes con lo que no hace falta calcularlo.

Otro ejemplo

qownnotes-media-cy7296

En este caso binomial con valores discretos, vemos que es imposible coger el 5%, lo más cercano sería 0.0352, no podremos elegir ni Zscore ni Ttest, por lo que tendremos que hechar mano del PValue

qownnotes-media-iN7296 qownnotes-media-wd7296

P-VALUE

Hemos visto cómo calcular comparar las media de dos grupos y decir cuál es su diferencia en un 95% de los casos con un inervalo de confianza, y hemos visto como contrastar los resultados de dos grupos, comporbando si el valor de uno estaba o no dentro del 5% de posibilidades del otro para poder descartar una hipótesis frente a la otra. Estos dos conceptos los unimos para porporcionar en un contraste de hipótesis el p-value, con el que compararmos nuestra hipótesis con la posibilidad de que al azar obtengamos ese mismo valor, proporcionando un intervalo de confianza.

Nos sirve para dictaminar que el resultado obtenido es estadísticamnete significativo.

qownnotes-media-AeDUXI

Es la probabilidad de obtener un resultado al menos tan extremo como el que realmente se ha obtenido (valor del estadístico calculado), suponiendo que la hipótesis nula es cierta. EL pvalue por tanto es una medida de significación estadística (un resultado o efecto es estadísticamente significativo cuando es improbable que haya sido debido al azar).

Una característica de p es que su valor depende del tamaño de la muestra a partir de la cual se determina: a mayor tamaño muestral mayor es la probabilidad de encontrar diferencias significativas entre los grupos que se comparan. De ahí la importancia de que en todo estudio deba realizarse el cálculo del tamaño muestral necesario en base a las caractrísticas de dicho estudio con anterioridad a la recogida y el análisis de los datos.

qownnotes-media-WX4036

Todos los libros de estadística tienen una tabla detrás con los pvalues asociado a cada uno de los valores de Zscore de distintas distribuciones, en este caso distribución normal

qownnotes-media-Tt4036

En primer paso por tanto para probar nuestra hipótesis será crear la H_o para tener una base sobre la que medir una alternativa a nuestra H_a con los datos obervados, calcular un test estadístico desde los datos dados, y finalmente comparar este test estadístico con la distribución hipotética. Esta comparación nos dice cómo de alejado el valor del test está de la hipótesis alternativa. Pvalue es la probabilidad sobre la H_0 de obtener evidencia tanto o más alejada de test estadístico obtenido de los datos observados en la dirección de la hipótesis alternativa.

Por tanto si pvalue es pequeño, o H_o es cierto y los datos observados tienen eventos raros o H_o es falso. Pvalye es el valor más pequeño de alfa en el cual rechazaríamos la H_0

P-value permite elegir un porcentaje de error Typo I , comparando ese porcenaje de error si P-value es menor rechazar la hipotesis nula. Hay muchos mal entendidos en su interpretación, es por eso que tiene muchos defensores y muchos detractores.

Vamos a calcular, cómo de inusal son nuestros dato estimado asumiento que Ho es cierto

qownnotes-media-PC7296

Para calcular pvalue en una distribución nomal, existe la función pnorm() con los distintos parámetros de la distribución, uno de ellos, lower.trail = TRUE/FALSE le decimos si queremos calcularlo por un sólo extremo o por los dos extremos.En el ejemplo queríamos ver el valor de alfa = 2, si sólo cogemos un extremo tendremos alfa >2

qownnotes-media-H11100

Si la ditribución es binomial pbinom.

qownnotes-media-wSCjdp

qownnotes-media-s11100 qownnotes-media-n11100

Si la distribución es de poisson ppois()

qownnotes-media-HOxTYi

qownnotes-media-u11100

Proporcionando el valor de pvalue en lugar de un alfa (probabilidad de cometer error tipo I) permites a los revisores de tu trabajo realizar hipótesis con el alpha que ellos elijan. La norma general será que si pvalue es menor que el alfa, rechazas H_o.

Si la hipótesis abarca ambos extremos, lo que se suele hacer es doblar el pvalue.

qownnotes-media-CPmAdQ

qownnotes-media-ptQNCu

qownnotes-media-dQbvin

qownnotes-media-tgvovm

qownnotes-media-krqBtz

qownnotes-media-QHoCgi

qownnotes-media-sNVxow

qownnotes-media-rMQxye

qownnotes-media-MKVoiM

qownnotes-media-kCHvKO

TEST

qownnotes-media-xOUDhE

qownnotes-media-uOtGOY

qownnotes-media-nLCkjG

qownnotes-media-SaFnex

qownnotes-media-ALBoXw

qownnotes-media-cPFSzZ

qownnotes-media-hbyEwA

PODER DE CONTRASTE

Poder es la probabilidad de rechazar la hipótesis nula cuando es falsa. El error tipo II es llamada Beta, es la proabilidad de no rechazar H_o cuando en realidad es falsa, por tanto el poder = 1-beta.

Se usa por ejemplo, cuando quieres determinar si el tamaño de la muestra es sufientemente grande para arrojar resultados significativos que no se consideren por azar.

En el ejemplo de eventos nocturnos que causan apnea del sueño, la zona azul que queda a la derecha de la línea vertical es el poder de contraste

qownnotes-media-C16176

El poder crecerá cuanto más grande sea la muestra y contra más alejado esté el velor que intentamos inferir del calculado por H_0, en el dibujo de arriba, cuando mayor sea la media de la curva azu, más se aleja del alfa que es la recta y negral mayor área caerá sobre ella a partir de la línea vertical que representa el 5% de H_0 (alfa o error tipo I). En este gráfico vemos cómo según aumenta el valor de media, el poder aumenta, y de manera mucho más incremental si aumentamos el número de muestras

qownnotes-media-sv4036

Cuanto más pequeá sea la desviación, más delgada será la campana, por lo que las curvas se solapan menos y su poder será mayor. osea poder de contraste decrecerá según aumente la desviación estándar ya que significa más variabilidad en los datos y más se solapen H_0 y H_A

qownnotes-media-Y16176

qownnotes-media-zvFTEF

Podemos calcular el poder de dos maneras, o bien manual o bien con la funcion t.test.

qownnotes-media-izEfWo

Para calcularlo de manera manual, usaremos la función pnorm():

z <- qnorm(.95) pnorm(30+z, mean = 32, lower.tail=FALSE)

Esto nos devolverá 0.64, que signiica que con un 64% podemos rechazar H_o cuando es falsa.

Tenemos la función power.t.test que nos realiza esta tarea para distribuciones Tstudent

qownnotes-media-fJ7448

En este caso vemos que tanto Delta (diferencia entre medias) como sd, si se cambian proporcionalmente, el cambio no afecta al poder. Este ration (Delta/sd) se llama effect size.

Se puede utilizar la misma fórmula para saber el número de individuos mínimos necesitamos para cumplir con un power concreto, esto es muy útil para saber el tamaño mínimo de la muestra

qownnotes-media-in7448

EJERCICIOS

Type swirl() when you are ready to begin.

Attaching package: ‘swirl’

The following objects are masked from ‘package:shinyjs’:

info, reset

Warning message: package ‘swirl’ was built under R version 3.3.3

swirl()

Welcome to swirl! Please sign in. If you’ve been here before, use the same name as you did then. If you are new, call yourself something unique.

What shall I call you? miguel

Would you like to continue with one of these lessons?

1: Statistical Inference Hypothesis Testing 2: Statistical Inference Probability2 3: Statistical Inference T Confidence Intervals 4: Statistical Inference Variance 5: No. Let me start something new.

Selection: 5

Please choose a course, or type 0 to exit swirl.

1: Exploratory Data Analysis 2: Getting and Cleaning Data 3: R Programming 4: Statistical Inference 5: Take me to the swirl course repository!

Selection: 4

Please choose a lesson, or type 0 to return to course menu.

1: Introduction 2: Probability1 3: Probability2 4: ConditionalProbability 5: Expectations
6: Variance 7: CommonDistros 8: Asymptotics 9: T Confidence Intervals 10: Hypothesis Testing
11: P Values 12: Power 13: Multiple Testing 14: Resampling

Selection: 12

Attempting to load lesson dependencies…

Package ‘reshape2’ loaded correctly!

Package ‘ggplot2’ loaded correctly!

| | 0%

Power. (Slides for this and other Data Science courses may be found at github https://github.com/DataScienceSpecialization/courses/. If you care to use them, they must be downloaded as a zip file and viewed locally. This lesson corresponds to 06_Statistical_Inference/11_Power.)

…

|= | 1%

In this lesson, as the name suggests, we’ll discuss POWER, which is the probability of rejecting the null hypothesis when it is false, which is good and proper.

…

|=== | 2%

Hence you want more POWER.

…

|==== | 3%

Power comes into play when you’re designing an experiment, and in particular, if you’re trying to determine if a null result (failing to reject a null hypothesis) is meaningful. For instance, you might have to determine if your sample size was big enough to yield a meaningful, rather than random, result.

…

|====== | 4%

Power gives you the opportunity to detect if your ALTERNATIVE hypothesis is true.

…

|======= | 5%

Do you recall the definition of a Type II error? Remember, errors are bad.

1: Rejecting a true null hypothesis 2: Miscalculating a t score 3: Misspelling the word hypothesis 4: Accepting a false null hypothesis

Selection: 4

All that practice is paying off!

|========= | 7%

Beta is the probability of a Type II error, accepting a false null hypothesis; the complement of this is obviously (1 - beta) which represents the probability of rejecting a false null hypothesis. This is good and this is POWER!

…

|========== | 8%

Recall our previous example involving the Respiratory Distress Index and sleep disturbances. Our null hypothesis H_0 was that mu = 30 and our alternative hypothesis H_a was that mu > 30.

…

|============ | 9%

Which of the following expressions represents our test statistic under this null hypothesis? Here X’ represents the sample mean, s is the sample std deviation, and n is the sample size. Assume X’ follows a t distribution.

1: (X’-30)/(s^2/n) 2: 30/(s/sqrt(n)) 3: (X’-30)/(s/sqrt(n)) 4: X’/(s^2/n)

Selection: 3

Nice work!

|============= | 10%

In the expression for the test statistic (X’-30)/(s/sqrt(n)) what does (s/sqrt(n)) represent?

1: a standard measure 2: a standard error 3: a standard sample 4: a standard variance 5: a standard bearer

Selection: 2

That’s the answer I was looking for.

|=============== | 11%

Suppose we’re testing a null hypothesis H_0 with an alpha level of .05. Since H_a proposes that mu > 30 (the mean hypothesized by H_0), power is the probability that the true mean mu is greater than the (1-alpha) quantile or qnorm(.95). For simplicity, assume we’re working with normal distributions of which we know the variances.

…

|================ | 12%

Here’s the picture we’ve used a lot in these lessons. As you know, the shaded portion represents 5% of the area under the curve. If a test statistic fell in this shaded portion we would reject H_0 because the sample mean is too far from the mean (center) of the distribution hypothesized by H_0. Instead we would favor H_a, that mu > 30. This happens with probability .05.

…

|================== | 13%

You might well ask, “What does this have to do with POWER?” Good question. We’ll look at some pictures to show you.

…

|=================== | 14%

First we have to emphasize a key point. The two hypotheses, H_0 and H_a, actually represent two distributions since they’re talking about means or centers of distributions. H_0 says that the mean is mu_0 (30 in our example) and H_a says that the mean is mu_a.

…

|===================== | 15%

We’re assuming normality and equal variance, say sigma^2/n, for both hypotheses, so under H_0, X’~ N(mu_0, sigma^2/n) and under H_a, X’~ N(mu_a, sigma^2/n).

…

|====================== | 16%

Here’s a picture with the two distributions. We’ve drawn a vertical line at our favorite spot, at the 95th percentile of the red distribution. To the right of the line lies 5% of the red distribution.

…

|======================== | 17%

Quick quiz! Which distribution represents H_0?

1: the blue 2: the red

Selection: 1

Nice try, but that’s not exactly what I was hoping for. Try again.

The two distributions have the same spread (variance). They differ in their means (centers). Which has a mean equal to that hypothesized by H_0?

1: the blue 2: the red

Selection: 2

Keep working like that and you’ll get there!

|========================= | 18%

Which distribution represents H_a?

1: the red 2: the blue

Selection: 2

Excellent work!

|=========================== | 20%

From the picture, what is the mean proposed by H_a?

1: 28 2: 32 3: 36 4: 30

Selection: 2

Nice work!

|============================ | 21%

See how much of the blue distribution lies to the right of that big vertical line?

…

|============================== | 22%

That, my friend, is POWER!

…

|=============================== | 23%

It’s the area under the blue curve (H_a) to the right of the vertical line.

…

|================================= | 24%

Note that the placement of the vertical line depends on the null distribution. Here’s another picture with fatter distributions. The vertical line is still at the 95th percentile of the null (red) distribution and 5% of the distribution still lies to its right. The line is calibrated to mu_0 and the variance.

…

|================================== | 25%

Back to our original picture.

…

|==================================== | 26%

We’ve shamelessly stolen plotting code from the slides so you can see H_a in action. Let’s look at pictures before we delve into numbers. We’ve fixed mu_0 at 30, sigma (standard deviation) at 4 and n (sample size) at 16. The function myplot just needs an alternative mean, mu_a, as argument. Run myplot now with an argument of 34 to see what it does.

myplot(34)

All that practice is paying off!

|===================================== | 27%

The distribution represented by H_a moved to the right, so almost all (100%) of the blue curve is to the right of the vertical line, indicating that with mu_a=34, the test is more powerful, i.e., there’s a higher probability that it’s correct to reject the null hypothesis since it appears false. Now try myplot with an argument of 33.3.

myplot(33.3)

You got it right!

|======================================= | 28%

This isn’t as powerful as the test with mu_a=34 but it makes a pretty picture. Now try myplot with an argument of 30.

myplot(30)

That’s the answer I was looking for.

|======================================== | 29%

Uh Oh! Did the red curve disappear? No. it’s just under the blue curve. The power now, the area under the blue curve to the right of the line, is exactly 5% or alpha!

…

|========================================== | 30%

So what did we learn?

…

|=========================================== | 32%

First, power is a function that depends on a specific value of an alternative mean, mu_a, which is any value greater than mu_0, the mean hypothesized by H_0. (Recall that H_a specified mu>30.)

…

|============================================= | 33%

Second, if mu_a is much bigger than mu_0=30 then the power (probability) is bigger than if mu_a is close to 30. As mu_a approaches 30, the mean under H_0, the power approaches alpha.

…

|============================================== | 34%

Just for fun try myplot with an argument of 28.

myplot(28)

You nailed it! Good job!

|================================================ | 35%

We see that the blue curve has moved to the left of the red, so the area under it, to the right of the line, is less than the 5% under the red curve. This then is even less powerful and contradicts H_a so it’s not worth looking at.

…

|================================================= | 36%

Here’s a picture of the power curves for different sample sizes. Again, this uses code “borrowed” from the slides. The alternative means, the mu_a’s, are plotted along the horizontal axis and power along the vertical.

…

|=================================================== | 37%

What does the graph show us about mu_a?

1: power is independent of mu_a 2: as it gets bigger, it gets less powerful 3: as it gets bigger, it gets more powerful

Selection: 3

All that hard work is paying off!

|==================================================== | 38%

What does the graph show us about sample size?

1: as it gets bigger, it gets more powerful 2: as it gets bigger, it gets less powerful 3: power is independent of sample size

Selection: 1

Excellent work!

|====================================================== | 39%

Now back to numbers. Our test for determining rejection of H_0 involved comparing a test statistic, namely Z=(X’-30)/(sigma/sqrt(n)), against some quantile, say Z_95, which depended on our level size alpha (.05 in this case). H_a proposed that mu > mu_0, so we tested if Z>Z_95. This is equivalent to X’ > Z_95 * (sigma/sqrt(n)) + 30, right?

…

|======================================================= | 40%

Recall that nifty R function pnorm, which gives us the probability that a value drawn from a normal distribution is greater or less than/equal to a specified quantile argument depending on the flag lower.tail. The function also takes a mean and standard deviation as arguments.

…

|========================================================= | 41%

Suppose we call pnorm with the quantile 30 + Z_95 * (sigma/sqrt(n)) and specify mu_a as our mean argument. This would return a probability which we can interpret as POWER. Why?

…

|========================================================== | 42%

Recall our picture of two distributions. 30 + Z_95 * (sigma/sqrt(n)) represents the point at which our vertical line falls. It’s the point on the null distribution at the (1-alpha) level.

…

|============================================================ | 43%

Study this picture. Calling pnorm with 30 + Z_95 * (sigma/sqrt(n)) as the quantile and mu_a, say 32, as the mean and lower.tail=FALSE does what?

1: returns the area under the red curve to the left of the line 2: returns the area under the red curve to the right of the line 3: returns the area under the blue curve to the left of the line 4: returns the area under the blue curve to the right of the line

Selection: 4

Your dedication is inspiring!

|============================================================= | 45%

Let’s try some examples now. Before we do, what do we know pnorm will return if we specify a quantile less than the mean?

1: an answer dependent on beta 2: an answer less than .50 3: an answer dependent on alpha 4: an answer greater than 1

Selection: 2

Excellent work!

|=============================================================== | 46%

First, define a variable z as qnorm(.95)

z <- qnorm(.95)

You are amazing!

|================================================================ | 47%

Run pnorm now with the quantile 30+z, mean=30, and lower.tail=FALSE. We’ve specified sigma and n so that the standard deviation of the sample mean is 1.

z <- pnorm(30+z,mean=30,lower.tail = FALSE)

Not quite, but you’re learning! Try again. Or, type info() for more options.

Type pnorm(30+z,mean=30,lower.tail=FALSE) at the command prompt.

pnorm(30+z,mean=30,lower.tail = FALSE) [1] 0.05

Your dedication is inspiring!

|================================================================== | 48%

That’s not surprising, is it? With the mean set to mu_0 the two distributions, null and alternative, are the same and power=alpha. Now run pnorm now with the quantile 30+z, mean=32, and lower.tail=FALSE.

pnorm(30+z,mean=32,lower.tail = FALSE) [1] 0.63876

You are really on a roll!

|=================================================================== | 49%

See how this is much more powerful? 64% as opposed to 5%. When the sample mean is quite different from (many standard errors greater than) the mean hypothesized by the null hypothesis, the probability of rejecting H_0 when it is false is much higher. That is power!

…

|==================================================================== | 50%

Let’s look again at the portly distributions.

…

|====================================================================== | 51%

With this standard deviation=2 (fatter distribution) will power be greater or less than with the standard deviation=1?

1: less than 2: the same 3: greater

Selection: 1

Keep up the great work!

|======================================================================= | 52%

To see this, run pnorm now with the quantile 30+z, mean=32 and sd=1. Don’t forget to set lower.tail=FALSE so you get the right tail.

pnorm(30+z,mean=32,sd=1,lower.tail = FALSE) [1] 0.63876

Excellent job!

|========================================================================= | 53%

Now run pnorm now with the quantile 30+z*2, mean=32 and sd=2. Don’t forget to set lower.tail=FALSE so you get the right tail.

pnorm(30+z,mean=32,sd=2,lower.tail = FALSE) [1] 0.5704709

Not exactly. Give it another go. Or, type info() for more options.

Type pnorm(30+z*2,mean=32,sd=2,lower.tail=FALSE) at the command prompt.

pnorm(30+z*2,mean=32,sd=2,lower.tail = FALSE) [1] 0.259511

Excellent job!

|========================================================================== | 54%

See the power drain from 64% to 26% ? Let’s review some basic facts about power. We saw before in our pictures that the power of the test depends on mu_a. When H_a specifies that mu > mu_0, then as mu_a grows and exceeds mu_0 increasingly, what happens to power?

1: it doesn’t change 2: it decreases 3: it increases

Selection: 3

You are quite good my friend!

|============================================================================ | 55%

Here’s another question. Recall our power curves from before.

…

|============================================================================= | 57%

As the sample size increases, what happens to power?

1: it decreases 2: it increases 3: it doesn’t change

Selection: 2

You nailed it! Good job!

|=============================================================================== | 58%

Here’s another one. More power curves.

…

|================================================================================ | 59%

As variance increases, what happens to power?

1: it doesn’t change 2: it increases 3: it decreases

Selection: 3

All that hard work is paying off!

|================================================================================== | 60%

Here’s another one. And even more power curves.

…

|=================================================================================== | 61%

As alpha increases, what happens to power?

1: it decreases 2: it doesn’t change 3: it increases

Selection: 3

Keep up the great work!

|===================================================================================== | 62%

If H_a proposed that mu != mu_0 we would calculate the one sided power using alpha / 2 in the direction of mu_a (either less than or greater than mu_0). (This is only approximately right, it excludes the probability of getting a large test statistic in the opposite direction of the truth.

…

|====================================================================================== | 63%

Since power goes up as alpha gets larger would the power of a one-sided test be greater or less than the power of the associated two sided test?

1: less than 2: greater 3: they’re the same

Selection: 3

Not quite! Try again.

The quantity alpha is bigger than alpha/2 so it’s got more power.

1: less than 2: greater 3: they’re the same

Selection: 2

All that hard work is paying off!

|======================================================================================== | 64%

Finally, if H_a specified that mu < mu_0 could we still do the same kind of power calculations?

1: No 2: Yes

Selection: 2

Nice work!

|========================================================================================= | 65%

Suppose H_a says that mu > mu_0. Then power = 1 - beta = Prob ( X’ > mu_0 + z_(1-alpha) * sigma/sqrt(n)) assuming that X’~ N(mu_a,sigma^2/n). Which quantities do we know in this statement, given the context of the problem? Let’s work through this.

…

|=========================================================================================== | 66%

What does the null hypothesis H_0 tell us that the population mean equals?

1: alpha 2: mu_0 3: beta 4: mu_a

Selection: 2

All that hard work is paying off!

|============================================================================================ | 67%

After the null mean mu_0 is proposed what does the designer of the hypothesis test specify in order to reject or fail-to-reject H_0? In other words, what is the level size of the test?

1: mu_0 2: alpha 3: mu_a 4: beta

Selection: 2

You nailed it! Good job!

|============================================================================================== | 68%

So we know that the quantities mu_0 and alpha are specified by the test designer. In the statement 1 - beta = Prob( X’ > mu_0 + z_(1-alpha) * sigma/sqrt(n)) given mu_a > mu_0, mu_0 and alpha are specified, and X’ depends on the data. The other four quantities, (beta, sigma, n, and mu_a), are all unknown.

…

|=============================================================================================== | 70%

It should be obvious that specifying any three of these unknowns will allow us to solve for the missing fourth. Usually, you only try to solve for power (1-beta) or the sample size n.

…

|================================================================================================= | 71%

An interesting point is that power doesn’t need mu_a, sigma and n individually. Instead only sqrt(n)*(mu_a - mu_0) /sigma is needed. The quantity (mu_a - mu_0) / sigma is called the EFFECT SIZE. This is the difference in the means in standard deviation units. It is unit free so it can be interpreted in different settings.

…

|================================================================================================== | 72%

We’ll work through some examples of this now. However, instead of assuming that we’re working with normal distributions let’s work with t distributions. Remember, they’re pretty close to normal with large enough sample sizes.

…

|==================================================================================================== | 73%

Power is still a probability, namely P( (X’ - mu_0)/(S /sqrt(n)) > t_(1-alpha, n-1) given H_a that mu > mu_a ). Notice we use the t quantile instead of the z. Also, since the proposed distribution is not centered at mu_0, we have to use the non-central t distribution.

…

|===================================================================================================== | 74%

R comes to the rescue again with the function power.t.test. We can omit one of the arguments and the function solves for it. Let’s first use it to solve for power.

…

|======================================================================================================= | 75%

We’ll run it three times with the same values for n (16) and alpha (.05) but different delta and standard deviation values. We’ll show that if delta (difference in means) divided by the standard deviation is the same, the power returned will also be the same. In other words, the effect size is constant for all three of our tests.

…

|======================================================================================================== | 76%

We’ll specify a positive delta; this tells power.t.test that H_a proposes that mu > mu_0 and so we’ll need a one-sided test. First run power.t.test(n = 16, delta = 2 / 4, sd=1, type = “one.sample”, alt = “one.sided”)$power .

power.t.test(n = 16, delta = 2 / 4, sd=1, type = “one.sample”,alt = “one.sided”)$power [1] 0.6040329

You are doing so well!

|========================================================================================================== | 77%

Now change delta to 2 and sd to 4. Keep everything else the same.

power.t.test(n = 16, delta = 2, sd=4, type = “one.sample”,alt = “one.sided”)$power [1] 0.6040329

Excellent job!

|=========================================================================================================== | 78%

Same answer, right? Now change delta to 100 and sd to 200. Keep everything else the same.

power.t.test(n = 16, delta = 100, sd=200, type = “one.sample”,alt = “one.sided”)$power [1] 0.6040329

Nice work!

|============================================================================================================= | 79%

So keeping the effect size (the ratio delta/sd) constant preserved the power. Let’s try a similar experiment except now we’ll specify a power we want and solve for the sample size n.

…

|============================================================================================================== | 80%

First run power.t.test(power = .8, delta = 2 / 4, sd=1, type = “one.sample”, alt = “one.sided”)$n .

power.t.test(power = .8, delta = 2 / 4, sd=1, type = “one.sample”, alt = “one.sided”)$n [1] 26.13751

Nice work!

|================================================================================================================ | 82%

Now change delta to 2 and sd to 4. Keep everything else the same.

power.t.test(power = .8, delta = 2, sd=4, type = “one.sample”, alt = “one.sided”)$n [1] 26.13751

Keep up the great work!

|================================================================================================================= | 83%

Same answer, right? Now change delta to 100 and sd to 200. Keep everything else the same.

power.t.test(power = .8, delta = 100, sd=200, type = “one.sample”, alt = “one.sided”)$n [1] 26.13751

Excellent job!

|=================================================================================================================== | 84%

Now use power.t.test to find delta for a power=.8 and n=26 and sd=1

power.t.test(power = .8, n = 26, sd=1, type = “one.sample”, alt = “one.sided”)

 One-sample t test power calculation

          n = 26
      delta = 0.5013986
         sd = 1
  sig.level = 0.05
      power = 0.8
alternative = one.sided

You’re close…I can feel it! Try it again. Or, type info() for more options.

Type power.t.test(power = .8, n=26, sd=1, type = “one.sample”, alt = “one.sided”)$delta.

power.t.test(power = .8, n = 26, sd=1, type = “one.sample”, alt = “one.sided”)$delta [1] 0.5013986

Great job!

|==================================================================================================================== | 85%

Not a surprising result, is it? It told you before that with an effect size of .5 and power .8, you need a sample size a little more than 26. Now run it with n=27.

power.t.test(power = .8, n = 27, sd=1, type = “one.sample”, alt = “one.sided”)$delta [1] 0.4914855

You’re the best!

|====================================================================================================================== | 86%

What do you think will happen if you doubled sd to 2 and ran the same test?

1: delta will halve 2: delta will double 3: delta won’t change

Selection: 1

Keep trying!

Since you’re doubling the denominator (sd) you have to double the numerator (delta) in order to keep the effect size constant.

1: delta won’t change 2: delta will halve 3: delta will double

Selection: 1

Not quite! Try again.

Since you’re doubling the denominator (sd) you have to double the numerator (delta) in order to keep the effect size constant.

1: delta will halve 2: delta won’t change 3: delta will double

Selection: 3

Great job!

|======================================================================================================================= | 87%

Now for a quick review. We call this the power.u.test since it comes after the power.t.test. LOL.

…

|========================================================================================================================= | 88%

1. The level of a test is specified by what?

1: gamma 2: None of the others 3: alpha 4: beta 5: delta

Selection: 5

One more time. You can do it!

The level refers to the probability of rejecting the null hypothesis when it is true. This is the first Greek letter.

1: gamma 2: delta 3: None of the others 4: alpha 5: beta

Selection: 4

Perseverance, that’s the answer.

|========================================================================================================================== | 89%

2. What is a Type II error?

1: accepting a true hypothesis 2: rejecting a false hypothesis 3: rejecting a true hypothesis 4: accepting a false hypothesis

Selection: 4

All that hard work is paying off!

|============================================================================================================================ | 90%

3. What is power?

1: delta 2: thrilling 3: None of the others 4: alpha 5: beta 6: gamma

Selection: 3

Perseverance, that’s the answer.

|============================================================================================================================= | 91%

4. You’re a perfectionist designing an experiment and you want both alpha and beta to be small. Can they both be 0 for this single test?

1: No 2: Yes

Selection: 1

Excellent job!

|=============================================================================================================================== | 92%

5. Suppose H_0 proposes mu = mu_0 and H_a proposes that mu < mu_0. You’ll test a series of mu_a with power != alpha. Which of the following is NOT true?

1: huh? 2: mu_a-mu_0 < 0 3: mu_a-mu_0=0 4: mu_0-mu_a > 0

Selection: 3

Keep up the great work!

|================================================================================================================================ | 93%

6. Suppose H_0 proposes mu = mu_0 and H_a proposes that mu < mu_0. Which of the following is true?

1: mu_0=mu_a maximizes the power 2: the smaller mu_0-mu_a the more powerful the test 3: the smaller mu_a-mu_0 the more powerful the test

Selection: 2

That’s not exactly what I’m looking for. Try again.

Here mu_a < mu_0 and the smaller mu_a-mu_0 is, the easier it is to discriminate between mu_a and mu_0.

1: the smaller mu_0-mu_a the more powerful the test 2: mu_0=mu_a maximizes the power 3: the smaller mu_a-mu_0 the more powerful the test

Selection: 3

You are doing so well!

|================================================================================================================================== | 95%

7. Which expression represents the size effect?

1: (mu_a - mu_0) / sqrt(sigma) 2: (mu_a - mu_0) / sqrt(n) 3: (mu_a - mu_0) / n 4: (mu_a - mu_0) / sigma

Selection: 4

You are doing so well!

|=================================================================================================================================== | 96%

8. True or False? More power is better than less power.

1: False 2: True

Selection: 2

That’s correct!

|===================================================================================================================================== | 97%

9. True or False? A larger beta (call it beta_max) is more powerful than a smaller beta.

1: True 2: False

Selection: 2

You’re the best!

|====================================================================================================================================== | 98%

10. True or False? The larger the sample size the less powerful the test.

1: False 2: True

Selection: 1

That’s correct!

|======================================================================================================================================== | 99%

Congrats! You finished this powerful lesson. We hope you feel emPOWERED.

…

|=========================================================================================================================================| 100%

Would you like to receive credit for completing this course on Coursera.org?

1: No 2: Yes

COMPARACIONES MÚLTIPLES (MULTIPLE TESTING)

Es un todo un subcampo detro de la inferencia estadística. Trata de reducir y ajustar los errores tipo I y II para evitar falsos positivos (rechazar H_o cuando no se debe) y falsos negativos (rechazar H_a cuando no se debe).

Imaginemos que queremos probar una hipótesis sobre múltiples grupos. Lo que haremos será comparar las medias obtenidas de todos esos grupos.

Un ejemplo claro que podemos tener en esta situación es que, pongamos que hacemos 10.000 test, si calculamos un PValue de manera que P< alfa estaremos controlando los falsos positivos en un nivel de media alfa. Suponemos P< 0.05 significante, por lo que 10000*0.05 = 500 posibles falsos positivos, que parece un número muy alto. Cómo evitamos esto?

Vamos a jugar con ALFA equilibrando la reducción del número de falsos positivos con que tengamos una ALFA muy baja y no podamos rechazar nunca Ho

qownnotes-media-y17684

FWER Para controlar el falso positivo: Family wise error rate (FWER) que es la probabilidad de encontrar al menos 1 falso positivo. Se puede intentar controlar este FWER con la corrección de Bonferroni. Lo bueno de este método es que es muy rápido, pero baja tanto ALFA que nos será muy dificil rechazar Ho

qownnotes-media-X17684

FDR

Otro método para controlar los falsos positivos es False discovery rate (FDR) most popular in genomics E(V/R) = Ratio de falsos positivos (número de falsos positivos divido número de positivos). En este caso se aplica la corrección BH (Benjamini-Hochberg). Es fácil de calcular y aunque permite más falsos positivos es menos conservador y permite rechazar más facilemnte Ho

qownnotes-media-Z17684

En esta gráfica, se ve que el número de resultados significantes sin corrección son 4, con la corrección HC baja a 3 y con la de bonferoni a 2

qownnotes-media-F17684

Adjusted Pvalues

Hay otra tercera aproximación que es en lugar de buscar un umbral a ALFA ajustar pvalues de las muestras

qownnotes-media-Y17684

qownnotes-media-fMlJWx

qownnotes-media-KNlCiF

qownnotes-media-WpNVqa

qownnotes-media-cGLCMd

qownnotes-media-ubeOBL

qownnotes-media-RrJxGj

En R hay una función llamada p.adjust que hace esto, teniendo un vector de mil pvalues:

qownnotes-media-d14368 qownnotes-media-c14368 qownnotes-media-r14368 qownnotes-media-w14368 qownnotes-media-Z14368 qownnotes-media-n14368

EJERCICIOS

Given that data is valuable and we’d like to get the most out of it, we might use it to test several hypotheses. If we have an alpha level of .05 and we test 20 hypotheses, then on average, we expect one error, just by chance.

…

|===== | 5%

Another potential problem is that after running several tests, only the lowest p-value might be reported OR all p-values under some threshold might be considered significant. Undoubtedly, some of these would be false.

…

|======= | 7%

Luckily, we have clever ways of minimizing errors in this situation. That’s what we’ll address. We’ll define specific error measures and then statistical ways of correcting or limiting them.

…

|========= | 8%

Multiple testing is particularly relevant now in this age of BIG data. Statisticians are tasked with questions such as “Which variables matter among the thousands measured?” and “How do you relate unrelated information?”

…

|========== | 10%

Since multiple testing addresses compensating for errors let’s review what we know about them. A Type I error is

1: rejecting a false hypothesis 2: failing to reject a true hypothesis 3: rejecting a true hypothesis 4: failing to reject a false hypothesis

Selection: 4

Nice try, but that’s not exactly what I was hoping for. Try again.

Eliminate the two choices that are not errors. A Type I error involves rejection.

1: rejecting a true hypothesis 2: failing to reject a true hypothesis 3: rejecting a false hypothesis 4: failing to reject a false hypothesis

Selection: 3

Keep trying!

Eliminate the two choices that are not errors. A Type I error involves rejection.

1: rejecting a true hypothesis 2: failing to reject a false hypothesis 3: failing to reject a true hypothesis 4: rejecting a false hypothesis

Selection: 1

All that practice is paying off!

|============ | 11%

In an American court, an example of a Type I error is

1: convicting an innocent person 2: letting the indicted off on a technicality 3: acquitting a guilty person

Selection: 3

Keep trying!

In an American court, the null hypothesis is that the accused is innocent. If he is convicted when he really is innocent then the null hypothesis is rejected incorrectly.

1: letting the indicted off on a technicality 2: acquitting a guilty person 3: convicting an innocent person

Selection: 3

You’re the best!

|============== | 13%

A Type II error is

1: rejecting a false hypothesis 2: failing to reject a false hypothesis 3: failing to reject a true hypothesis 4: rejecting a true hypothesis

Selection: 1

Almost! Try again.

Eliminate the two choices that are not errors. A Type II error involves failing to reject.

1: failing to reject a false hypothesis 2: rejecting a true hypothesis 3: rejecting a false hypothesis 4: failing to reject a true hypothesis

Selection: 1

That’s correct!

|=============== | 15%

In an American court, an example of a Type II error is

1: letting the indicted off on a technicality 2: convicting an innocent person 3: acquitting a guilty person

Selection: 3

Great job!

|================= | 16%

Good. Let’s continue reviewing. The null hypothesis

1: tells us the origins of the number 0 2: represents the status_quo and is assumed true 3: is a big nothing that statisticians like to gossip about 4: is never true

Selection: 2

That’s the answer I was looking for.

|=================== | 18%

The p-value is “the probability under the null hypothesis of obtaining evidence as or more extreme than your test statistic (obtained from your observed data) in the direction of the alternative hypothesis." Of course p-values are related to significance or alpha levels, which are set before the test is conducted (often at 0.05).

…

|==================== | 20%

If a p-value is found to be less than alpha (say 0.05), then the test result is considered statistically significant, i.e., surprising and unusual, and the null hypothesis (the status quo) is ?

1: revised 2: renamed the aleph null hypothesis 3: rejected 4: accepted

Selection: 3

You got it right!

|====================== | 21%

Now consider this chart copied from http://en.wikipedia.org/wiki/Familywise_error_rate. Suppose we’ve tested m null hypotheses, m_0 of which are actually true, and m-m_0 are actually false. Out of the m tests R have been declared significant, that is, the associated p-values were less than alpha, and m-R were nonsignificant, or boring results.

…

|======================== | 23%

Looking at the chart, which variables are known?

1: A,B,C 2: S,T,U,V 3: m_0, and m 4: m and R

Selection: 4

That’s a job well done!

|========================== | 25%

In testing the m_0 true null hypotheses, V results were declared significant, that is, these tests favored the alternative hypothesis. What type of error does this represent?

1: Type I 2: Type III 3: Type II 4: a serious one

Selection: 1

Excellent work!

|=========================== | 26%

Another name for a Type I error is False Positive, since it is falsely claiming a significant (positive) result.

…

|============================= | 28%

Of the m-m_0 false null hypotheses, T were declared nonsignificant. This means that these T null hypotheses were accepted (failed to be rejected). What type of error does this represent?

1: Type III 2: Type II 3: a serious one 4: Type I

Selection: 2

Perseverance, that’s the answer.

|=============================== | 30%

Another name for a Type II error is False Negative, since it is falsely claiming a nonsignificant (negative) result.

…

|================================ | 31%

A rose by any other name, right? Consider the fraction V/R.

…

|================================== | 33%

The observed R represents the number of test results declared significant. These are ‘discoveries’, something different from the status quo. V is the number of those falsely declared significant, so V/R is the ratio of FALSE discoveries. Since V is a random variable (i.e., unknown until we do an experiment) we call the expected value of the ratio, E(V/R), the False Discovery Rate (FDR).

…

|==================================== | 34%

A rose by any other name, right? How about the fraction V/m_0? From the chart, m_0 represents the number of true H_0’s and m_0 is unknown. V is the number of those falsely declared significant, so V/m_0 is the ratio of FALSE positives. Since V is a random variable (i.e., unknown until we do an experiment) we call the expected value of the ratio, E(V/m_0), the FALSE POSITIVE rate.

…

|====================================== | 36%

Another good name for the false positive rate would be

1: false alarm rate 2: the Type II rate 3: a thorn 4: a rose

Selection: 1

You are really on a roll!

|======================================= | 38%

The false positive rate would be closely related to

1: the Type I error rate 2: the Type II error rate 3: a thorny rose

Selection: 1

Perseverance, that’s the answer.

|========================================= | 39%

We call the probability of at least one false positive, Pr(V >= 1) the Family Wise Error Rate (FWER).

…

|=========================================== | 41%

So how do we control the False Positive Rate?

…

|============================================ | 43%

Suppose we’re really smart, calculate our p-values correctly, and declare all tests with p < alpha as significant. This means that our false positive rate is at most alpha, on average.

…

|============================================== | 44%

Suppose we perform 10,000 tests and alpha = .05. How many false positives do we expect on average?

1: 50000 2: 5000 3: 50 4: 500

Selection: 4

Keep working like that and you’ll get there!

|================================================ | 46%

You got it! 500 false positives seems like a lot. How do we avoid so many?

…

|================================================= | 48%

We can try to control the family-wise error rate (FWER), the probability of at least one false positive, with the Bonferroni correction, the oldest multiple testing correction.

…

|=================================================== | 49%

It’s very straightforward. We do m tests and want to control the FWER at level alpha so that Pr(V >= 1) < alpha. We simply reduce alpha dramatically. Set alpha_fwer to be alpha/m. We’ll only call a test result significant if its p-value < alpha_fwer.

…

|===================================================== | 51%

Sounds good, right? Easy to calculate. What would be a drawback with this method?

1: too many results will fail 2: too many results will pass 3: requires too much math

Selection: 1

That’s the answer I was looking for.

|======================================================= | 52%

Another way to limit the false positive rate is to control the false discovery rate (FDR). Recall this is E(V/R). This is the most popular correction when performing lots of tests. It’s used in lots of areas such as genomics, imaging, astronomy, and other signal-processing disciplines.

…

|======================================================== | 54%

Again, we’ll do m tests but now we’ll set the FDR, or E(V/R) at level alpha. We’ll calculate the p-values as usual and order them from smallest to largest, p_1, p_2,…p_m. We’ll call significant any result with p_i <= (alpha*i)/m. This is the Benjamini-Hochberg method (BH). A p-value is compared to a value that depends on its ranking.

…

|========================================================== | 56%

This is equivalent to finding the largest k such that p_k <= (k * alpha)/m, (for a given alpha) and then rejecting all the null hypotheses for i=1,…,k.

…

|============================================================ | 57%

Like the Bonferroni correction, this is easy to calculate and it’s much less conservative. It might let more false positives through and it may behave strangely if the tests aren’t independent.

…

|============================================================= | 59%

Now consider this chart copied from the slides. It shows the p-values for 10 tests performed at the alpha=.2 level and three cutoff lines. The p-values are shown in order from left to right along the x-axis. The red line is the threshold for No Corrections (p-values are compared to alpha=.2), the blue line is the Bonferroni threshold, alpha=.2/10 = .02, and the gray line shows the BH correction. Note that it is not horizontal but has a positive slope as we expect.

…

|=============================================================== | 61%

With no correction, how many results are declared significant?

1: 8 2: 6 3: 2 4: 4

Selection: 4

Excellent job!

|================================================================= | 62%

With the Bonferroni correction, how many tests are declared significant?

1: 6 2: 2 3: 4 4: 8

Selection: 2

All that practice is paying off!

|================================================================== | 64%

So the Bonferroni passed only half the results that the No Correction (comparing p-values to alpha) method passed. Now look at the BH correction. How many tests are significant with this scale?

1: 1 2: 7 3: 5 4: 3

Selection: 4

You’re the best!

|==================================================================== | 66%

So the BH correction which limits the FWER is between the No Correction and the Bonferroni. It’s more conservative (fewer significant results) than the No Correction but less conservative (more significant results) than the Bonferroni. Note that with this method the threshold is proportional to the ranking of the values so it slopes positively while the other two thresholds are flat.

…

|====================================================================== | 67%

Notice how both the Bonferroni and BH methods adjusted the threshold (alpha) level of rejecting the null hypotheses. Another equivalent corrective approach is to adjust the p-values, so they’re not classical p-values anymore, but they can be compared directly to the original alpha.

…

|======================================================================== | 69%

Suppose the p-values are p_1, … , p_m. With the Bonferroni method you would adjust these by setting p’_i = max(m * p_i, 1) for each p-value. Then if you call all p’_i < alpha significant you will control the FWER.

…

|========================================================================= | 70%

To demonstrate some of these concepts, we’ve created an array of p-values for you. It is 1000-long and the result of a linear regression performed on random normal x,y pairs so there is no true significant relationship between the x’s and y’s.

…

|=========================================================================== | 72%

Use the R command head to see the first few entries of the array pValues.

head(pValues) [1] 0.5334915 0.2765785 0.8380943 0.6721730 0.8122037 0.4078675

Great job!

|============================================================================= | 74%

Now count the number of entries in the array that are less than the value .05. Use the R command sum, and the appropriate Boolean expression.

sum() [1] 0

That’s not the answer I was looking for, but try again. Or, type info() for more options.

Type sum(pValues < 0.05) at the command prompt.

sum(pValues < 0.05) [1] 51

You are doing so well!

|============================================================================== | 75%

So we got around 50 false positives, just as we expected (.05*1000=50). The beauty of R is that it provides a lot of built-in statistical functionality. The function p.adjust is one example. The first argument is the array of pValues. Another argument is the method of adjustment. Once again, use the R function sum and a boolean expression using p.adjust with method=“bonferroni” to control the FWER.

p.adjust(pValue,method=“bonferroni”) Error in p.adjust(pValue, method = “bonferroni”) : object ‘pValue’ not found p.adjust(pValues,method=“bonferroni”)

sum(p.adjust(pValues,method=“bonferroni”) < 0.05) [1] 0

Your dedication is inspiring!

|================================================================================ | 77%

So the correction eliminated all the false positives that had passed the uncorrected alpha test. Repeat the same experiment, this time using the method “BH” to control the FDR.

sum(p.adjust(pValues,method=“BH”) < 0.05) [1] 0

Perseverance, that’s the answer.

|================================================================================== | 79%

So the BH method also eliminated all the false positives. Now we’ve generated another 1000-long array of p-values, this one called pValues2. In this data, the first half ( 500 x/y pairs) contains x and y values that are random and the second half contain x and y pairs that are related, so running a linear regression model on the 1000 pairs should find some significant (not random) relationship.

…

|==================================================================================== | 80%

We also created a 1000-long array of character strings, trueStatus. The first 500 entries are “zero” and the last are “not zero”. Use the R function tail to look at the end of trueStatus.

tail(trueStatus) [1] “not zero” “not zero” “not zero” “not zero” “not zero” “not zero”

Great job!

|===================================================================================== | 82%

Once again we can use R’s greatness to count and tabulate for us. We can call the R function table with two arguments, a boolean such as pValues2<.05, and the array trueStatus. The boolean obviously has two outcomes and each entry of trueStatus has one of two possible values. The function table aligns the two arguments and counts how many of each combination (TRUE,“zero”), (TRUE,“not zero”), (FALSE,“zero”), and (FALSE,“not zero”) appear. Try it now.

table(trueStatus) trueStatus not zero zero 500 500

That’s not the answer I was looking for, but try again. Or, type info() for more options.

Type table(pValues2 < 0.05, trueStatus) at the command prompt.

table(pValues2 < 0.05, trueStatus) trueStatus not zero zero FALSE 0 476 TRUE 500 24

You are doing so well!

|======================================================================================= | 84%

We see that without any correction all 500 of the truly significant (nonrandom) tests were correctly identified in the “not zero” column. In the zero column (the truly random tests), however, 24 results were flagged as significant.

…

|========================================================================================= | 85%

What is the percentage of false positives in this test?

4.8 [1] 4.8

Give it another try. Or, type info() for more options.

Divide 24 by 500 to get the percentage.

24/500 [1] 0.048

Your dedication is inspiring!

|========================================================================================== | 87%

Just as we expected - around 5% or .05*100.

…

|============================================================================================ | 89%

Now run the same table function, however, this time use the call to p.adjust with the “bonferroni” method in the boolean expression. This will control the FWER.

table(pValues2 < 0.05, trueStatus) trueStatus not zero zero FALSE 0 476 TRUE 500 24

Nice try, but that’s not exactly what I was hoping for. Try again. Or, type info() for more options.

Type table(p.adjust(pValues2,method=“bonferroni”) < 0.05, trueStatus) at the command prompt.

table(p.adjust(pValues2,method=“bonferroni”) < 0.05, trueStatus) trueStatus not zero zero FALSE 23 500 TRUE 477 0

You nailed it! Good job!

|============================================================================================== | 90%

Since the Bonferroni correction method is more conservative than just comparing p-values to alpha all the truly random tests are correctly identified in the zero column. In other words, we have no false positives. However, the threshold has been adjusted so much that 23 of the truly significant results have been misidentified in the not zero column.

…

|=============================================================================================== | 92%

Now run the same table function one final time. Use the call to p.adjust with “BH” method in the boolean expression. This will control the false discovery rate.

table(p.adjust(pValues2,method=“BH”) < 0.05, trueStatus) trueStatus not zero zero FALSE 0 487 TRUE 500 13

Great job!

|================================================================================================= | 93%

Again, the results are a compromise between the No Corrections and the Bonferroni. All the significant results were correctly identified in the “not zero” column but in the random (“zero”) column 13 results were incorrectly identified. These are the false positives. This is roughly half the number of errors in the other two runs.

…

|=================================================================================================== | 95%

Here’s a plot of the two sets of adjusted p-values, Bonferroni on the left and BH on the right. The x-axis indicates the original p-values. For the Bonferroni, (adjusting by multiplying by 1000, the number of tests), only a few of the adjusted values are below 1. For the BH, the adjusted values are slightly larger than the original values.

…

|===================================================================================================== | 97%

We’ll conclude by saying that multiple testing is an entire subfield of statistical inference. Usually a basic Bonferroni/BH correction is good enough to eliminate false positives, but if there is strong dependence between tests there may be problems. Another correction method to consider is “BY”.

…

|====================================================================================================== | 98%

Congrats! We hope you liked the multiple concepts and questions you saw in this lesson.

…

|========================================================================================================| 100%

RESAMPLING INFERENCE

Hasta ahora hemos visto métodos de inferencia estadística para obtener los parámetros de una población en base a su distrubución, pero ¿qué pasa si desconocemos su distribución?¿Cómo podemos obterner la media y varianza poblacional por ejemplo al tirar un dado si sabemos que está trucado? De lo que se trata entonces es de realizar inferencias sobre nuestros datos tomando muestras de nuestros propios datos, es decir haciendo múltiples selecciones aleatorias sobre los datos que hemos recogido, para poder obterner la media de sus medias.

BOOTSTRAPING

Es una técnica que permite simplificar la obtención de los intervalos de confianza en lugar de con métodos matemáticos, calcularlos a través de múltiples muestras aleatorias tomadas de las muestras sin sacar del espacio muestral en cada iteración la muesta obtenida (reemplazo).

Es una técnica no paramétrica y nos serviría para poder detallar en nuestro estudio un intervalo de confianza y el error estándar sobre una única medición tomada, iterando una y otra vez sobre esa muesta para obtener de manera aleatoria uno se sus valores y así simular la población.

qownnotes-media-JkhSIH

qownnotes-media-tanpiX

qownnotes-media-oTTWHa

qownnotes-media-HmdAlZ

qownnotes-media-MvzZKg

qownnotes-media-ugudYv

qownnotes-media-CGFVKT

qownnotes-media-zblFDp

PERMUTATION TEST

Esta es otra técnica de inferencia por muestreo, en la que en lugar de utilizar T-Student o Chicuadrado, como no sabemos la distribución de los datos comparamos los grupos haciendo permutaciones aleatorias entre ellos mismos y lo comparamos con los resultados más extremos que podemos obtener de su comparativa inicial.

A grandes rasgos, se restan las medias de ambos grupos y se comprueba con otro grupo de permutaciones realizadas cuántas restas de las medias de las permutaciones aleatorias supera la medida inicial, y así saber cuán extremo es la diferencia.

qownnotes-media-WVBDaB

qownnotes-media-lrLfbM

qownnotes-media-RRTRbe

qownnotes-media-BxnNue

qownnotes-media-txJbfA

EJERCICIOS

The bootstrap is a handy tool for making statistical inferences. It is used in constructing confidence intervals and calculating standard errors for statistics that might be difficult for some reason (e.g., lack of data or no closed form). Wikipedia tells us that bootstrapping is a technique which “allows estimation of the sampling distribution of almost any statistic using very simple methods." Simple is good, right?

…

|==== | 4%

The beauty of bootstrapping is that it avoids complicated mathematics and instead uses simulation and computation to infer distributional properties you might not otherwise be able to determine.

…

|====== | 6%

It’s relatively new, developed in 1979, by Bradley Efron, a Stanford statistician. The basic bootstrap principle uses OBSERVED data to construct an ESTIMATED population distribution using random sampling with replacement. From this distribution (constructed from the observed data) we can estimate the distribution of the statistic we’re interested in.

…

|======= | 7%

So, in bootstrapping the observed data substitutes for what?

1: a statistic 2: a hypothesis 3: a population 4: observations

Selection: 1

Not quite right, but keep trying.

We’ll used the observed data and sample it with replacement, just as we would do to find out information about some population.

1: a statistic 2: a population 3: observations 4: a hypothesis

Selection: 2

You nailed it! Good job!

|========= | 8%

So, in bootstrapping if the observed data is the population, what would the random samplings correspond to?

1: a statistic 2: a hypothesis 3: observations 4: a population

Selection: 3

You got it right!

|========== | 10%

In effect, the original observed sample substitutes for the population. Our samplings become observations from which we estimate a statistic and get an idea about its distribution. This lets us better understand the underlying population (from which we didn’t have enough data).

…

|============ | 11%

Here’s a critical point. In constructing the estimated distribution we sample the observed data WITH replacement. If the original sample is n long and we sampled n times without replacement what would we get?

1: an entirely new sample 2: a worse sample 3: the original sample permuted 4: a better sample

Selection: 3

Great job!

|============= | 12%

The motivating example from the slides involves computing the average of 50 rolls of a die. Of course we can do this theoretically when we know that the die is fair. Remember, E(x) = Sum(x*p(x)) for x=1,2,…6, and p(x)=1/6 for all values of x.

…

|============== | 14%

For the heck of it, compute the expected die roll for a fair die.

3 [1] 3

Give it another try. Or, type info() for more options.

Type sum(1:6)/6 at the command prompt.

sum(1:6)/6 [1] 3.5

Perseverance, that’s the answer.

|================ | 15%

Theoretically, the average is 3.5. Here, we’ve run code and plotted a histogram after we took 1000 such averages, each of 50 dice rolls. Note the unusual y-axis scale. We’re displaying this as a density function so the area of the salmon-colored region is theoretically 1. With this scale, though, all the heights of the bins actually add up to 5. So you have to multiply each height by .2 and add up all the results to get 1.

…

|================= | 17%

The point is, the empirical matches the theoretical. Yay! The highest bin is centered at 3.5 just as the math predicted. So what?

…

|=================== | 18%

What if some joker wanted you to run the same experiment with a die he gave you and he warned you that the dice was loaded? In other words, it wasn’t fair. It has some random distribution like this.

…

|==================== | 19%

The outcomes aren’t equally likely, are they? So when you do your 1000 runs of 50 rolls each, the density of the means looks different.

…

|====================== | 21%

We’ve done this for you and put the result in g2. Type print(g2) now to see the picture.

print(g2)

You’re the best!

|======================= | 22%

Picture’s a little different, right? Although this example is a bit contrived, it illustrates an important concept. We really want a distribution of means and we have only one set of observations. (In this case it was the empirical distribution associated with the unfair die - the big blue picture.) We used that one distribution, to “create” many (1000) distributions by sampling with replacement from the given one. We sampled 50000 times so we created 1000 distributions of 50 rolls each.

…

|========================= | 24%

We then calculated the mean of each of our created distributions and got a distribution of means. Sampling the one distribution many times gives us some variability in the resulting statistics we calculate. We can then calculate the standard error and confidence intervals associated with the statistic.

…

|========================== | 25%

Before we go on to more theory, here’s another example in which we’ll try to find a distribution of medians of a population. Do you recall what a median is?

1: the most frequent outcome 2: 50th percentile 3: a person who talks to spirits 4: a point halfway between rare and well-done

Selection: 1

Not quite right, but keep trying.

Half the outcomes are above it and half below it.

1: the most frequent outcome 2: a person who talks to spirits 3: 50th percentile 4: a point halfway between rare and well-done

Selection: 3

That’s the answer I was looking for.

|=========================== | 26%

Recall the father and son height data. Once again, we’ve loaded it for you. We’ve placed the height of the sons in the vector sh and the length of this vector is stored in the variable nh. Use the R command head to look at the first few entries of sh.

head(sh) [1] 59.77827 63.21404 63.34242 62.79238 64.28113 64.24221

That’s correct!

|============================= | 28%

Now look at nh to see how long sh is.

length(sh) [1] 1078

That’s not the answer I was looking for, but try again. Or, type info() for more options.

Type nh at the command prompt.

nh [1] 1078

Keep working like that and you’ll get there!

|============================== | 29%

Now we’ll create 1000 distributions of the same length as the original sh. We’ll do this by sampling sh with replacement 1000*nh times and store the results in an array with 1000 rows, each with nh entries. Then we’ll take the median of each row and plot the result.

…

|================================ | 31%

Note that every time we draw from the empirical distribution sh, each of its nh data points is equally likely to be pulled, therefore the probability of drawing any one is 1/nh. The 1000 samples we create will vary from the original.

…

|================================= | 32%

Here’s the resulting density curve. This estimates the distribution of medians. The thick vertical line shows where the median of the original, observed data sh lies.

…

|=================================== | 33%

We stored the 1000 medians of the resampled sets in the vector resampledMedians. Use the R function median to compute the median of numbers in this vector.

mean(nh) [1] 1078

Not quite! Try again. Or, type info() for more options.

Type median(resampledMedians) at the command prompt.

median(resampledMedians) [1] 68.61273

You got it!

|==================================== | 35%

Now compute the median of the original sample sh.

warning messages from top-level task callback ‘mini’ Warning message: In rm(list = c(vName), envir = globalenv()) : object ‘resamples’ not found > mean(sh) [1] 68.68407

Not quite! Try again. Or, type info() for more options.

Type median(sh) at the command prompt.

median(sh) [1] 68.61582

All that hard work is paying off!

|====================================== | 36%

Pretty close, right? Now back to theory. Suppose you have a statistic that estimates some population parameter, but you don’t know its sampling distribution. The bootstrap principle uses the distribution defined by the observed data to approximate the sampling distribution of that statistic.

…

|======================================= | 38%

The nice thing about bootstrapping is that you can always do it with simulation. The general procedure follows by first simulating B complete data sets from the observed data by sampling with replacement. Make sure B is large and that you’re sampling WITH replacement to create data sets the same size as the original.

…

|======================================== | 39%

This approximates drawing from the sampling distribution of that statistic, at least as well as the data approximates the true population distribution. By calculating the statistic for each simulated data set and using these simulated statistics we can either define a confidence interval (e.g. find the 2.5 and the 97.5 percentiles) or take the standard deviation to estimate a standard error of that statistic.

…

|========================================== | 40%

Notice that this process doesn’t use any fancy math or asymptotics. The only assumption behind it is that the observed sample is representative of the underlying population.

…

|=========================================== | 42%

We’ve created the vector fh for you which contains the fathers’ heights from the father son data we’ve been working with. It’s the same length as the sons’ data (1078) which is stored in nh. B, the number of bootstraps we want has been set to 1000. We’ll do an example now in small steps.

…

|============================================= | 43%

Our one sample of observed data is in the vector fh. Use the R function sample to sample fh nh*B times. Set the argument replace to TRUE. Put the result in the variable sam.

sam <-sample(fh,nh*B,replace = TRUE)

All that hard work is paying off!

|============================================== | 44%

Now form sam into a matrix with B rows and nh columns. Use the R function matrix and put the result in resam

resam <- matrix

Almost! Try again. Or, type info() for more options.

Type resam <- matrix(sam,B,nh) at the command prompt.

resam <- matrix(sam,B,nh)

You are quite good my friend!

|================================================ | 46%

Now use the R function apply to take the median (third argument) of each row of resam (first argument). Put the result in meds. The second argument, the number 1, specifies that the application of the function is to the rows of the first argument.

meds <- apply(resam,1,median)

Excellent work!

|================================================= | 47%

| Now look at the difference between the median of fh and the median of meds. > median(meds)-median(fh + ) [1] -0.004105

You are amazing!

|=================================================== | 49%

Pretty close, right? Now use the R function sd to estimate the standard error of the vector meds.

sd(meds) [1] 0.1041539

You are quite good my friend!

|==================================================== | 50%

We previously did this same process for the sons’ data and stored the resampled medians in the 1000-long vector resampledMedians. Find the standard error of resampledMedians.

sd(resampledMedians) [1] 0.08293445

You got it right!

|===================================================== | 51%

Now we’ll find a 95% confidence interval for the sons’ data with the R function quantile. The first argument is the vector of resampledMedians and the second is the expression c(.025,.975). Do this now.

quantile(resampledMedians,c(.025,.975)) 2.5% 97.5% 68.43733 68.80718

You are amazing!

|======================================================= | 53%

Pretty close quantiles, right? Now do the same thing for the fathers’ data. Recall that it’s stored in the vector meds.

quantile(meds,c(.025,.975)) 2.5% 97.5% 67.54999 67.94110

You’re the best!

|======================================================== | 54%

Another pair of close quantiles, but notice that these quantiles of the fathers’ medians differ from those of the sons.

…

|========================================================== | 56%

Bootstrapping is a very diverse and complicated topic and we just skimmed the surface here. The technique we showed you is nonparametric, that is, it’s not based on any parameterized family of probability distributions. We used only one set of observations that we assumed to be representative of the population.

…

|=========================================================== | 57%

Finally, the confidence intervals we calculated might not perform very well because of biases but the R package bootstrap provides an easy fix for this problem.

…

|============================================================= | 58%

Now, to permutation testing, another handy tool used in group comparisons. As bootstrapping did, permutation testing samples a single dataset a zillion times and calculates a statistic based on these samplings.

…

|============================================================== | 60%

Permutation testing, however, is based on the idea of exchangability of group labels. It measures whether or not outcomes are independent of group identity. Our zillion samples simply permute group labels associated with outcomes. We’ll see an example of this.

…

|================================================================ | 61%

Here’s a picture from the dataset InsectSprays which contains counts of the number of bugs killed by six different sprays.

…

|================================================================= | 62%

We’ll use permutation testing to compare Spray B with Spray C.

…

|================================================================== | 64%

Use the R command dim to find the dimensions of InsectSprays.

dim(InsectSprays) [1] 72 2

That’s a job well done!

|==================================================================== | 65%

Now use the R command names to find what the two columns of InsectSprays contain.

names(InsectSprays) [1] “count” “spray”

You got it!

|===================================================================== | 67%

We’ll use permutation testing to compare Spray B with Spray C. We subsetted data for these two sprays into a data frame subdata. Moreover, the two data frames Bdata and Cdata contain the data for their respective sprays.

…

|======================================================================= | 68%

Now use the R command range on Bdata$count to find the minimum and maximum counts for Spray B.

range(Bdata$count) [1] 7 21

Excellent work!

|======================================================================== | 69%

The picture makes more sense now, right? Now do the same for Spray C. Its data is in Cdata.

range(Cdata$count) [1] 0 7

All that practice is paying off!

|========================================================================== | 71%

From the ranges (as well as the picture), the sprays look a lot different. We’ll test the (obviously false) null hypothesis that their means are the same.

…

|=========================================================================== | 72%

To make the analysis easier we’ve defined two arrays for you, one holding the counts for sprays B and C. It’s call BCcounts. Look at it now.

BCcounts [1] 11 17 21 11 16 14 17 17 19 21 7 13 0 1 7 2 3 1 2 1 3 0 1 4

Excellent job!

|============================================================================= | 74%

The second array we’ve defined holds the spray identification and it’s called group. These two arrays line up with each other, that is, the first 12 entries of counts are associated with spray B and the last 12 with spray C. Look at group now.

group [1] “B” “B” “B” “B” “B” “B” “B” “B” “B” “B” “B” “B” “C” “C” “C” “C” “C” “C” “C” “C” “C” “C” “C” “C”

You got it!

|============================================================================== | 75%

We’ve also defined for you a one-line function testStat which takes two parameters, an array of counts and an array of associated identifiers. It assumes all the counts come from group B or group C. It subtracts the mean of the counts from group C from the mean of the counts of group B. Type testStat with no parentheses and no arguments to see how it’s defined.

testStat() Error in mean(w[g == “B”]) : argument “w” is missing, with no default testStat function(w, g) mean(w[g == “B”]) - mean(w[g == “C”]) <environment: 0x0000000021419b48>

Excellent job!

|=============================================================================== | 76%

Now set a variable obs by invoking testStat with the arguments BCcounts and group and assigning the result to obs.

obs <- testStat((BCcounts,group)) Error: unexpected ‘,’ in “obs <- testStat((BCcounts,” obs <- testStat(BCcounts,group)

Great job!

|================================================================================= | 78%

Take a peek at obs now.

obs [1] 13.25

Nice work!

|================================================================================== | 79%

Pretty big difference, right? You can check this by using mean on Bdata$count and on Cdata$count and subtracting the latter from the former. Equivalently, you can just apply mean to Bdata$count-Cdata$count. Do either one now.

apply(Bdata$count-Cdata$count,1,mean) Error in apply(Bdata$count - Cdata$count, 1, mean) : dim(X) must have a positive length apply(Bdata$count-Cdata$count,mean) Error in match.fun(FUN) : argument “FUN” is missing, with no default mean(Bdata$count-Cdata$count) [1] 13.25

You are really on a roll!

|==================================================================================== | 81%

So, mean(Bdata$count)-mean(Cdata$count) equals mean(Bdata$count-Cdata$count) because ?

1: mean is linear 2: the data is special 3: mathemagic

Selection: 1

Nice work!

|===================================================================================== | 82%

Now this is where the permutation testing starts to involve resampling. We’re going to test whether or not the particular group association of the counts affects the difference of the means.

…

|======================================================================================= | 83%

We’ll keep the same array of counts, just randomly relabel them, by permuting the group array. R makes this process very easy. Calling the function sample (which we’ve used several times in this lesson) with one argument, an array name, will simply permute the elements of that array.

…

|======================================================================================== | 85%

Call sample now on the array group to see what happens.

sample(gruop) Error in sample(gruop) : object ‘gruop’ not found sample(group) [1] “C” “C” “B” “C” “C” “B” “C” “B” “C” “B” “B” “B” “C” “B” “B” “C” “C” “B” “C” “C” “B” “C” “B” “B”

You nailed it! Good job!

|========================================================================================== | 86%

The labels are all mixed up now. We’ll do this permuting of labels and then we’ll recalculate the difference of the means of the two “new” (really newly labelled) groups.

…

|=========================================================================================== | 88%

We’ll relabel and calculate the difference of means 10000 times and store the differences (of means) in the array perms. Here’s what the code looks like perms <- sapply(1 : 10000, function(i) testStat(BCcounts, sample(group))). Try it now.

perms <- sapply(1 : 10000, function(i) testStat(BCcounts,sample(group)))

Great job!

|============================================================================================ | 89%

We can take the mean of the virtual array of the boolean expression perms > obs. Do this now.

mean(pems > pems) Error in mean(pems > pems) : object ‘pems’ not found mean(perms > obs) [1] 0

Nice work!

|============================================================================================== | 90%

So on average 0 of the permutations had a difference greater than the observed. That means we would reject the null hypothesis that the means of the two sprays were equal.

…

|=============================================================================================== | 92%

Here’s a histogram of the difference of the means. Looks pretty normal, right? We can see that the distribution runs roughly between -10 and +10 and it’s centered around 0. The vertical line shows where the observed difference of means was and we see that it’s pretty far away from the distribution of the resampled permutations. This means that group identification did matter and sprays B and C were quite different.

…

|================================================================================================= | 93%

Here’s the picture of the InsectSprays again. Suppose we run the same experiment, this time comparing sprays D and E, which look more alike. We’ve redefined testStat to look at these sprays and subtract the mean of spray E from the mean of spray D.

…

|================================================================================================== | 94%

We’ve also stored off the D and E data in DEcounts and the group labels in group. Run testStat now with DEcounts and group.

testStat(DEcounts,group) [1] 1.416667

That’s correct!

|==================================================================================================== | 96%

We’ve stored off this value, 1.416667, in the variable obs for you. Now run the permutation command, with DEcounts. Here it is, perms <- sapply(1 : 10000, function(i) testStat(DEcounts, sample(group)))

perms <- sapply(1 : 10000, function(i) testStat(DEcounts, sample(group)))

Great job!

|===================================================================================================== | 97%

Finally, we can plot the histogram of the distribution of the difference of the means. We see that with these sprays the observed difference of means (the vertical line) is closer to the mean of the permuted labels. This indicates that sprays D and E are quite similar and we fail to reject the null hypothesis that the means were equal.

…

|======================================================================================================= | 99%

Congrats! We hope you weren’t bugged too much by this lesson and feel like you’ve pulled yourself up by your bootstraps.

TEST

qownnotes-media-MNbivq

qownnotes-media-fFZicc

qownnotes-media-VMctaN

qownnotes-media-DwoZEa

qownnotes-media-oXujCs

qownnotes-media-pPrnhV

qownnotes-media-GcRabm

qownnotes-media-ETEsxE

qownnotes-media-RnXkAC